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ABSTRACT 


This report summarizes the research performed by North Carolina State University and 
NASA Ames Research Center under Cooperative Agreement NCA2-7 19, "Numerical 
Simulation of Supersonic and Hypersonic Inlet Flow Fields". Four distinct rotated upwind 
schemes were developed and investigated to determine accuracy and practicality. The scheme 
found to have the best combination of attributes, including reduction to grid alignment with no 
rotation, was the cell centered non-orthogonal (CCNO) scheme. In 2-D, the CCNO scheme 
dramatically improved accuracy when used with first order flux interpolation. CCNO also 
improved rotation when flux interpolation was extended to second order. In 3-D, 
improvements were less dramatic in all cases, with second order flux interpolation showing the 
least improvement over grid aligned upwinding. The reduction in improvement is attributed to 
uncertainly in determining optimum rotation angle and difficulty in performing accurate and 
efficiently interpolation angle in 3-D. The CCNO rotational technique will prove very useful 
for increasing accuracy when second order interpolation is not appropriate and will materially 
improve inlet flow solutions. 


INTRODUCTION 

The development programs of the High Speed Civil Transport, National Aerospace 
Plane, and next generation fighter/attack aircraft have identified the need for greater 
understanding and predictive capability of the complex fluid dynamic processes occurring in 
high speed inlet configurations. The flow fields in the inlet are replete with intricate fluid 
phenomena such as crossing shock waves, shock-wave/turbulent-boundary-layer interactions, 
merging comer flow boundary layers, and vortical flow. Because of limitations in wind tunnel 
capabilities such as limitations of scale and in non-obtrusive measuring techniques, much of the 
understanding of inlet flowfields is a result of numerical simulation. These computer solutions 
rely heavily on state-of-the-art numerical techniques. 

A commonly and confidently used numerical procedure is the upwind scheme which 
has been a major breakthrough in the modeling of fluid flow in the transonic through 
hypersonic flow regimes. The upwind schemes have provided the ability to capture flow 
discontinuities within a few grid points without tuning the artificial dissipation. However, it is 
well known that excess dissipation is generated by upwind schemes when the captured 
discontinuities are oblique to the computational grid. The excess dissipation smears the 
captured shock waves, thereby reducing the accuracy of the numerical predictions. High speed 
inlet configurations are especially prone to this effect because the grid topologies are naturally 
h-type configurations with oblique crossing shock waves, reflecting shock waves, and shock- 
wave/boundary-layer interactions. Thus it is the objective of this research to develop a 
numerical model that improves the prediction and resolution of these dominant flow features in 
regions where they are captured oblique to the computational mesh and thereby improve the 
accuracy of the global inlet solution. 

A promising approach to reducing the aforementioned smearing effect is employing a 
rotated upwind scheme. The rotated upwind procedure is one of a class of newly developing 
numerical procedures that attempt to introduce multi-dimensional dependency into the 
numerical integration. The basic idea of the rotated upwind scheme is to dynamically align the 
upwind difference stencil in a direction based on the developing flow field features. This 
procedure is an alternative to traditional schemes where the dissipation model is affixed by the 
computational mesh. The realignment of the upwind operator enhances the ability of the 
computational model to predict more accurately the local flow physics. The rotation is 
especially beneficial in regions where a dominant flow feature, such as a shock wave, exists 
and is oblique to the computational grid. As mentioned previously, these features commonly 
occur in high speed inlet simulations. 



A focused investigation was conducted under this cooperative agreement to determine 
how the rotation could best be performed and to ascertain how rotation improves the 
computational results. Rather than repeat published explanations and data, a review of the 
techniques developed and results obtained is presented below and two AIAA papers which 
include detailed explanations and results are included as appendices. These papers mclude 
brief reviews of previous rotation and multidimensional upwind research, including references. 

RESEARCH RESULTS 

Several rotated upwind schemes have been developed for the Euler equations in two 
dimensions, all of which show significantly improved shock capturing ability over grid aligned 
schemes to first-order accuracy when the shock wave is oblique to the grid. It is initially 
unclear which of the schemes offer the best promise for further development. Therefore an 
initial study is performed where the essential differences between the previous schemes are 
distilled and then used as building blocks for competing baseline algorithms. This step of the 
research resulted in four, first-order accurate, rotated upwind schemes that are presented in 
AIAA Paper 94-0079 "Rotated Upwind Strategies for Solution of the Euler Equations", a copy 
of which is attached. The four baseline schemes are the result of a two parameter survey of 
algorithm characteristics. The parameters are the position in the cell at which the rotation takes 
place (cell-center vs. cell-edge) and the topological space in which the rotation takes place 
(physical vs. computational). It is found that the cell-center rotation strategies are more robust 
and accurate than the cell-edge schemes because of a requisite averaging procedure in the cell- 
center schemes. The averaging procedure acts as an inherent smoothing agent and eliminates 
odd-even decoupling which is found to be characteristic of the cell-edge schemes. It is also 
shown that rotating in computational space vs. a rotation in physical space is simpler to 
implement and has the advantage that it will revert to a grid aligned formulation. It is also 
shown that the first-order accurate rotated scheme produces results typical of the second-order 
accurate grid aligned scheme for the solution of a two-dimensional reflecting shock wave duct 
configuration. 

The comparative study noted above revealed that a baseline first-order accurate rotated 
upwind scheme where the rotation in computational space is performed at each cell-center and 
designated as CCNO (Cell Center rotation, Non-Orthogonal fluxes) has the best overall 
characteristics. The scheme was then extended to second-order accuracy and viscous terms 
were included to model the Navier- Stokes equations. Furthermore, rotated boundary 
conditions were also developed. The extended algorithm was applied to various shock 
reflection and shock-boundary layer interaction flowfields. The results are contained in AIAA 
Paper 94-2291, "Rotated Upwind Algorithms for Solution of Two- and Three-Dimensional 
Euler and Navier-Stokes Equations" which is attached. The first test case, that of an inviscid 
shock wave reflection demonstrates that the CCNO algorithm yields a more accurate result than 
the grid aligned scheme to both first- and second-order accuracy. Furthermore, the improved 
accuracy is maintained as the grid density is increased. This result is an improvement over 
previous rotated schemes which show only marginal improvement to second-order accuracy. 
The next test case is a shock wave impinging on a turbulent boundary layer. To first-order 
accuracy, the CCNO scheme is shown to produce results in better agreement with the 
experimental data as compared to the grid aligned scheme. To second-order accuracy, although 
the inviscid region of the flowfield is qualitatively improved with the CCNO scheme, no 
improvement in the wall pressure and skin friction distributions are realized. This result is not 
unexpected, as the mathematical character of the governing equations changes from essentially 
hyperbolic outside the viscous layer to essentially parabolic near the surface. Therefore despite 
the improved prediction in the inviscid region, the wall properties are dominated by the viscous 
terms prediction which are inherently elliptical in nature and thereby do not benefit from 
upwinding. Although this is a disappointing result in terms of wall predictions, it is 



encouraging in terms of interior flow measurements such as mass flow rates and thrust 
performance which are integrated quantities across the entire domain. 

The two-dimensional CCNO algorithm was extended to three dimensions and is also 
described in the previously cited paper. The key development made during the extension of the 
algorithm was the derivation of two sequences of rotation that align a coordinate axis in any 
given preferred direction. Furthermore, the sequences of coordinate rotation are designed to 
lake advantage of directional symmetry such that the number of possible orientations of the 
rotated coordinate system with respect to the original coordinate system is minimized. 
Consequently, the logic required in the interpolation and projection procedures is reduced. 
The result is an automatic method of aligning efficiently a coordinate axis in any given or 
computed preferred direction. 

The three-dimensional algorithm is shown to be sufficiently robust to compute complex 
flowfields commonly found in supersonic through hypersonic inlet configurations. An inviscid 
three-dimensional shock surface reflection test case shows that the CCNO algorithm improves 
three-dimensional shock wave capturing to both first- and second-order accuracy. This is a 
significant result in light of other recent multi-dimensional algorithms that show only marginal 
or no improvement to three-dimensional inviscid calculations. However, it is also shown that 
the accuracy improvements in the three-dimensional calculation are not as great as what was 
seen in two dimensions. It is believed that this is a result of uncertainties in selecting a true 
dominant direction which becomes more problematic in three dimensions as compared to two 
dimensions. Viscous calculations of an intersecting-wedge^comer-flow geometry and a 
generic hypersonic inlet configuration (the inlet results are contained in Kontinos, D. A. 
"Rotated Upwind Algorithms for Solution of the Two- and Three-Dimensional Euler and 
Navier-Stokes Equations," Ph.D. Dissertation, Department of Mechanical and Aerospace 
Engineering, North Carolina State University, Raleigh, NC, 1994) show that as in the two- 
dimensional viscous calculations, the inviscid portion of the flowfield, e.g., the shock waves, 
are qualitatively improved using the CCNO algorithm. However to first-order accuracy, the 
predictions of pressure and heat transfer on the wall are not improved using the rotated 
algorithm. To second-order accuracy, there is some improvement in the wall property 
properties in regions of a dominant cross flow feature. For instance, in the intersecting wedge 
case better agreement with the experimental data is achieved with the CCNO algorithm in the 
cross-flow reattachment region of the multiple shock wave structure in the corner of the 
geometry. Similar improvement is seen in the generic inlet configuration in the region where 
the boundary layer is rolled into a comer vortex. 

The accuracy improvements gained by the CCNO algorithm come at a considerable 
cost Both the two- and three-dimensional algorithms are shown to consume 2-4 times more 
computer time than their grid aligned counteiparts. The increase in computer usage is caused 
by an increase in the work required per iteration and a reduced stability limit which causes an 
increase in the number of iterations for convergence. However, regions of potential algorithm 
improvement are highlighted that may reduce the work load, such as a more accurate temporal 
linearization and a simplification of the interpolation procedure. It is believed that the CCNO 
scheme may have application for inviscid flow solutions but the scheme is immature for 
viscous flowfields. In viscous flowfields, selection of a dominant direction becomes uncertain 
in subsonic regions such as in boundary layers, especially in three dimensions. Moreover, the 
fine grid spacing typical of viscous flow field meshes restricts and allowable time step and 
increases the amount of computer usage. As mentioned previously, a more accurate temporal 
linearization may overcome this difficulty. Finally, current research into truly multi- 
dimensional Riemann solvers may reveal an optimum upwind direction in the boundary layer 
thereby refining the CCNO schemes effectiveness. 



CONCLUSIONS 


The development and evaluation of rotated upwind procedures revealed that both 
accuracy and practicality must be considered in the choice of such schemes. The rotation 
scheme that was chosen, CCNO, was the best compromise between accuracy and ease of 
implementation. This scheme was demonstrated to improve results in essentially mviscid flow 
reeimes in both 2-D and 3-D. The improvement was most dramatic when a first order accurate 
rotation was used, producing results that compared well with second-order accurate grid 
aligned solutions. The second order CCNO scheme also gave improved results over second 
order grid aligned results. However, the improvement was not as dramatic in 3-D and was 
achieved at considerable computational cost. Further research is indicated to determine how 
best to find the optimum upwinding direction in 3-D and how to obtain efficiently interpolated 
flux quantities in the rotated directions. The lack of improvement due to rotation in the viscous 
layer is forecast by the mathematical character of the governing equations in this region. 

RECOMMENDATIONS 


CCNO may be used to increase accuracy to match second order grid aligned results 
with only first order flux interpolation. This may be particularly useful when the gnd is 
distorted to conform to geometry or flow features. 

Further research should be conducted to optimize rotation direction and to improve 
higher order flux interpolation, particularly in 3-D. 
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Abstract 

Four finite volume, fully conservative, rotated upwind 
strategies for solution of the two-dimensional Euler 
equations are compared. The four strategies are based on the 
combinations of two options which are a cell-edge vs. a 
cell-center rotation, and a rotation in physical space vs. 
computational space. The four schemes are implemented 
with maximum commonalty for direct comparison. The 
solutions are relaxed to steady state by the LU-SGS scheme. 
Solutions of a Mach 2 channel flow problem are presented 
and compared in terms of accuracy and robustness. It is 
shown that the cell-edge strategies create unacceptable 
oscillations in the solution while the cell-center strategies 
contain inherent smoothing that allow for accurate solutions 
with good convergence properties. The first-order rotated 
upwind results are seen to be comparable to standard second- 
order upwind results. 

1 fl Introduction 

The development of upwind schemes has been a major 
breakthrough in the modeling of fluid flow in the transonic 
through hypersonic regimes. The upwind schemes have 
provided the ability to capture flow discontinuities within a 
few grid points without tuning the artificial dissipation. 
However, it is well known that excess dissipation is 
generated when the captured discontinuities are oblique to 
the grid. 

One such upwind algorithm. Roe's 1 scheme, has been 
used quite successfully by a variety of researchers. The Roe 
c^hpmet models flow discontinuities as a series of linearized 
waves. In one dimension, these linearized waves model the 
actual wave propagation quite accurately. However, in 
multidimensions, the one dimensional Roe scheme 
misinterprets the multidimensional flow field and 
incorrectly models local wave propagation. This manifests 
itself as excess dissipation and results in the smearing of 
flow discontinuities. An excellent discussion of this effect 
is presented in Ref. [2], 
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In an effort to improve upwind schemes in 
pi yiHHimpnsinnal flow, much recent effort has been focused 
on Hi«fgipatinn models that incorporate multidimensional 
effects. A comprehensive review by van Leer 3 gives a good 
overview of such methods. One of the approaches is the 
use of a rotated upwind scheme. The basic idea is to orient 
the upwind solver in a preferred direction (or 'rotate the 
integration stencil to this preferred direction) such that the 
one dimensional operator can more accurately model the 
multidimensional flow. Several rotated upwind schemes 
have been developed for the Euler equations in two 
dimensions^* 3, all of which show significant 
im provement over grid aligned schemes to first-order 
accuracy and modest improvement to second-order accuracy. 
The main drawback of these algorithms has been the 
increase in computer work as compared to the grid aligned 
algorithms. Moreover, to date a three dimensional 
algorithm has not been developed and it is unclear whether a 
rotated formulati on will show the same improvements in 
3D as it does in 2D. 

The objective of this study is to compare four two- 
dimensional rotated upwind strategies in terms of accuracy 
and robustness so that a method can be selected for 
development in three dimensions. Each of the strategies 
contain elements of previously developed algorithms, 
although several aspects are new. They are implemented 
with as much commonalty as possible to promote a good 
comparison. 

This study is organized as follows; first Section 2 
presents the governing equations. Section 3 presents a 
general overview of each of the four strategies. The details 
of the implementation of the strategies are presented in 
Section 4. Section 5 will present results of a simple 
rhannid flow problem and Section 6 will conclude. 

2 0 Thfl Flllef Equations 

For this study, the governing equations are the two- 
dimensional Euler equations in integral conservation law 
form. The set in index notation is as follows: 

— jpdV+ Jp« i n 1 dS = 0 , mass 

m Vol Surf 

— Jp Uj dV+ J (pu,U; + pS^rti dS = 0 , mom 

* Vol Smf 



d_ 

dt 


J pE,dV+ | (E, + p)Mi«i dS = Q . 

Vol Surf 


energy 


In the finite volume framework, the conserved variables are 
taken as averages over the cell volume. Thus, the conserved 
quantities can be taken out of the volume integrals in the 
previous equations. Moreover, the surface integrals become 
summations over the cell sides. Expanding out the 
momentum equation into its vector components, the 
equations can be written as, 


method wherein a preferred upwind direction is selected at 
each cell face and two orthogonal fluxes are computed; an 
upwind flux in the preferred direction, and a central 
difference flux normal to the preferred direction. This is 
shown schematically in Fig. 1. Later versions of the 
algorithm apply Roe's scheme in both directions^. The 
values of the dependent variables are interpolated between 
cell centers through a linear-quadratic interpolating function. 
The rotated fluxes are then projected onto the grid face 
direction to ensure conservation. Levy shows improved 
shock capturing over grid aligned schemes. 


dU 

dt 


-J-Yf.s . 

Vol * 


sides 


where the conserved variables are, 

U = [p pu pv E f ] r , 


( 2 . 1 ) 


( 22 ) 


and the fluxes are given by, 

F=[pV puV + pi pvV + pj (E, + p)v] , (2.3) 

where the cell face normal with magnitude equal to the cell 
face length is denoted by n g and the velocity by V . The 
variables p , p, u, and v are the pressure, density, and x and y 
Parthian velocity components, respectively. Using the 
perfect gas relation, the total energy E, is given by, 


The «""" strategy is taken by Dadone and Grossman 6 . 
However, a different tactic is used to compute the numerical 
flux. An upwind flux (Roe's scheme) is computed in both 
the primary (preferred) and secondary (orthogonal) direction. 
Also a simplified interpolation stencil is incorporated that 
essentially 'grabs' the values of the characteristic variables at 
the nearest cell center to the rotated stencil. Computations 
using this yhpinp. also show improved results as compared 
to grid aligned schemes. Moreover, convergence rate and 
boundary condition improvements are shown as compared to 
previous rotated schemes. 

The cell-edge, orthogonal flux philosophy is 
incorporated in this study and is termed CEO. An upwind 
flux is computed in both the preferred and secondary 
direction at each cell edge. The interpolation of the 
primitive variables is linear between adjoining cells. 

3 2 Cell-Center Orthogona l Flux Strategy fCCO) 


E,= 


P 

(r-i) 


p(u 2 + v 2 ) 


(2.4) 


where y is the ratio of specific heats which has a value of 
1.4 for air. 


3 0 Rotate d T Tpwind Strategies 


This section presents an overview of the four strategies 
compared in this study. Recall that the main goal of a 
rotated upwind scheme is to align the one dimensional 
Riemann solver across flow discontinuities to reduce excess 
dissipation. In order to achieve this goal, one must 
compute the flux in a direction independent of the grid. 
This section will describe four strategies for such a flux 
computation. Each strategy will base the flux computation 
at the cell edge or cell center, and compute the preferred 
direction in physical space or computational space. Ibis 
section describes the strategies in broad terms while Section 
4 presents the details. 

3.1 Cell-Edge Orthogonal Flux Strategy (CEO) 

The first strategy is a distillation of the scheme 
developed by Levy, et. al. 5 . Levy's scheme is a cell-edge 


As an alternative to the previous cell-edge strategy, the 
concept of a cell-center rotation is now proposed (the notion 
of assigning a rotation angle to a cell center is recognized 
by Davis 4 although a method is not developed). This 
second strategy, termed cell-center orthogonal or CCO, is 
derived from the previous rotated finite difference scheme of 
Kontinos and McRae 7 . In this previous work, a rotation 
angle is selected in computational space thereby defining a 
rotated c oordinate system. The flux divergence is then 
computed in the rotated frame. The method shows 
improved results over grid aligned schemes for both in viscid 
and viscous flows. Unfortunately, the grid point residuals 
are computed independently of each other and the scheme is 
not conservative as originally implemented. 

In order to guarantee conservation, the scheme is 
converted from finite difference to finite volume to become 
the CCO strategy. In the CCO scheme, four fluxes are 
computed in two orthogonal directions originating at the 
cell center (recall in the cell-edge philosophy, two 
orthogonal fluxes are computed at each cell face and then 
projected onto that cell face). These four fluxes are then 
projected around the cell edges, details of which are 
presented in the next section. The CCO strategy is 
represented in Fig. 2 where a single cell is shown 
schematically along with the rotated fluxes. Because each 
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cell is rotated independently, non-unique fluxes at the cell 
edges are created. More precisely, since each cell edge is 
adjacent to two cell centers, two values of the flux are 
computed. For example the flux at face i+l/2 is computed 
from the rotation at cell i and i+i (see Fig. 2). To uniquely 
define the flux, the projected fluxes from each of the 
adjacent cell centers are averaged. 

A consequence of the CCO scheme is that for arbitrary 
grids, the scheme will not revert to a grid aligned 
formulation. From Fig. 2 one can see that a set of 
orthogonal axes will not align with a cell whose sides are 
not orthogonal. 

An important comparison between the CEO and CCO 
strategies can now be made. In two dimensions, the two 
schemes require about the same amount of arithmetic 
operations (work). Let W be the amount of work for n 
points. The work is mostly composed of the intetpolation 
and the flux computation, and is given by, 


Fig. 3a where a set of orthogonal axes are rotated through 
an angle 0 in computational space. Mapped into physical 
space, the fluxes are computed in non-orthogonal 
directions as shown in Fug. 3b. Thus this strategy is 
termed cell-center non-orthogonal or CCNO. These fluxes 
are projected around the cell as in the previous CCO 
strategy. However, in this case since the fluxes are not 
orthogonal, the projection must be performed by a 
coordinate transformation which is described in Section 4. 

3.4 Cell-Edge Non-Orthogonal Flux Strategy (CENQ) 

The ideas of the previous strategy are coupled with the 
cell-edge method to lead to the cell-edge non-orthogonal 
(CENO) strategy. As in the CCNO technique, orthogonal 
axes are rotated in computational space to a preferred 
direction with the cell edge acting as the origin. This 
results in two non-orthogonal fluxes being computed in 
physical space that are projected onto the grid face. The 
method is shown schematically in Fig. 4. 
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2 fluxes + 4inte ip^ 2 faces xnpo . nts 
J point 


face 

= 4n fluxes + 8n interp , 
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4 fluxes + 4 intern . 

xn points 

face J 


= 4n fluxes + 4n interp. 


However, an advantage of the cell-center scheme over the 
cell-edge scheme is predicted in three dimensions for an 
unrestricted rotation. In three dimensions, the cell-edge 
strategy requires three fluxes computed at each cell face. 
The cell-center strategy requires six fluxes at each cell 
center. The work is now given by. 


( 3 fluxes + 6 interp^ , faces . , 

Wrr = \ x3 xn points 

V face J point 

= 9 n fluxes + 18n interp , 


W, 


cc 


-(■ 


6 fluxes + 6 interp 


face 

= 6n fluxes + 6n interp 


x n points 


So it is seen that the cell-center strategy potentially reduces 
the work load by at least 1/3. 


4.0 Implementation 

This section provides the details of the implementation 
of each of the four strategies. They are implemented with 
as much commonalty as possible. First, all fluxes are 
computed using Roe’s scheme. Second, linear interpolation 
of the primitive variables between cell centers is used to 
obtain the state values for the flux computation. The 
interpolated values are used to compute both the inviscid 
flux and the upwind damping term (Riemann solution). It 
is possible to use higher order interpolations such as a one- 
dimensional quadratic function, or two-dimensional bilinear, 
reduced biquadratic, or biquadratic functions. But, most of 
these functions admit overshoots or undershoots in the 
interpolating function. Thus, the simplest interpolation is 
used to compare the four strategies. The order of accuracy is 
determined by the interpolation stencil which for this study 
is first-order accurate. Third, the rotation angles are 
computed in an analogous fashion between strategies. Each 
cell edge or center is assigned its own rotation angle which 
is not averaged with surrounding angles. Finally, each 
strategy is relaxed to steady state using the LU-SGS scheme 
of Yoon and Jameson 10 . 

Sections 4.2-4.S presents the details for the three steps 
peculiar to the rotated flux computation: the rotation angle 
calculation, the interpolation, and the flux projection. 


3.3 Cell-Center Non-Orthogonal Flux Strategy (CCNO) 4 . 1 Genera li zed Pro j ect i on 


The desire to have the CCO scheme revert to a grid 
aligned mode leads to the third method. This again is 
derived from the previous work of Kontinos and McRae 7 . 
In this method, the cell center rotation is performed in 
computational space where for rotation angles of n/i/2, the 
scheme automatically aligns the difference stencil in the grid 
contravariant directions. The rotated stencil is shown in 


In order to guarantee conservation, each of the four 
methods require the rotated fluxes to be projected back onto 
the grid faces. Thus, a generalized projection is presented. 
The development in this section is taken from Ref. [11]. 
Consider a field point as in Fig. 5. Five vectors are 
shown: h g is a vector in the direction of the cell face 
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normal with magnitude equal to the cell face length; n and 
n 2 are two arbitrary non-collinear vectors in the direction 
the flux is to be computed and are designated as the 
contravariant directions; and and /ij are the covanant 
base vectors. The covariant vectors are related to the 
contravariant vectors by. 


where. 


n 1 = ^n 2 xi), ='Jg(kxn 1 ) , ( 41 > 

-L = (s‘xn 2 ).Jk . (4.2) 


The flux vector expressed in the covariant vector base is 
given by. 


F = (F«n 1 )n 1 +(F*n 2 )ft 2 . 


The flux vector projected onto the cell face (see Eq. (2.1)) is 
given by, 


^ + ^ x i+ij+in + ^ x tj-u* + ) ' 

= *4 + ^yi+KJ+vj + + ^Ti+kj-in ) * 

and the gradient is, 

<4 - 8) 

A rotation direction can be computed directly from the 
gradient through the inverse tangent function. However, in 
order to eliminate oscillations of the rotation angle in the 
freestream, the following blending between the computed 
gradient direction 9^ and the grid face normal 9 g is 

incorporated, 

• (45) 


)«! *n f +(F*n 2 )« 2 *V 

(4.4) 

( -DM IN) 


- - - — * T71..M Atl 


c? 

*> 

II 

3 

(4.10) 


The details of the computation of the rotation angle, 
interpolation, and projection will be presented for an i+1/2 
face in the CEO strategy. An analogous computation is 
performed for the j+1/2 faces and will be omitted for 
brevity. 

The rotation angle is based on a flow field gradient. 
Let a flow field quantity be represented by q. The local 
gradient is computed as, 

4rigbt = +9i+l,/+l)> 

<Zieft = + 2«u + Qui* 1 ) . (4 - 5) 

Qap = T , (?I,;+1 "*■ Qi+\,j+\ )> 

9down — 1 9i+l,;-l)’ 


left. W 

J| = 2( 9up ” ^dow,, )' 

Then locally at face i+1/2 the metrics are given as, 

= tU2J ’ = ’ (4.1) 


Then the cosine of the rotation angle 9 is given as, 

cos(d) = <» cos(0 ? ) + (1 - tu) cos( 9 g ) . (4. 1 1) 

The quantity on acts as a blending coefficient controlled by 
the input parameter DMIN. 

The interpolation of the primitive variables is 
performed by mapping the local stencil into (%,V) space 
where -0.5££:S0.5 and -IStjSI. This is shown 
schematically in Fig. 6. The mapping is based on the 
averaged metrics of Eq. (4.7). Although holding the metrics 
of the transformation constant over the interpolation stencil 
is not strictly accurate in general, it is adequate for grids 
where the metrics do not change drastically. From Fig. 6 it 
is seen that four (|,rj) pairs are needed to compute the 
interpolation. A point on the interpolating ray is given by, 
£' = &> +Z x dx + J y dy, (4.12) 

V'=V o + rixdx+ri y dy. 

Setting the origin to (0,0) and using the rotation angle, Eq. 
(4.12) becomes, _ 

= § x cos(0) + sin(0), (4.13) 

q' = Tj x cos(9) + v y sm{9). 

The value of (£\tj’) must be clipped to lie in the 
interpolating region. One option would be to have (£rj) lie 
on the ellipse drawn in Fig. 6 and use a bilinear, or linear- 
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quadratic interpolating function over the local two- 
dimensional region. The option used in this study is to clip 
the interpolating ray to lie on the rectangular stentil of the 
surrounding cell centers and use linear interpolation. The 
clipping is given by. 


In this case, the Jacobian of the transformation is unity and 
the set of contravariant and covariant vectors are identical. 
Equation (4.4) becomes, 

F»n t =(F*n 1 )cos{0-0,)-(F*n2)sm(0-fl,).(4.22) 


SW = 


[l for2||is|n'| 

10 for 2 ||'| < 


5 = sign(l') 


i (1 STV) 111 l 


T/ = sign(j 7 ')| 1-SW + 


SW 


V ] 

2 (|£'| + l-Sw) 


(4.14) 

(4.15) 


The quantity SW is a switching function and is included in 
the denominators of Eq. (4.15) to avoid a divide by zero. 
The interpolant value denoted by q is then given by. 


q = c 0 + a tj (-c 0 + Crfj+iji + c 2 qiji), (4.16) 


where, 

q = ±+ 1 , c 2 = i - <* , c 0 = c x q M%i + c 2 q,j , (4.17) 


f 1 for tj i 0 
|-1 for tj < 0 


(4.18) 


jl = j + a. 


(4.19) 


This interpolation function is used for all four (£,7j) pairs of 
Fig. 6. The second (£tj) pair is given by, 

= -1 X sin(0) + f, cos(0) , (4.20) 

tj' = -7 j x sin(fl) + rj y cos( 0 ) , 


then the clipping function of Eq. (4.15) is used. The final 
two pairs are simply given by (< 5 , Jj ) 3 ^ 

(^ t 7 j ) 4 = (-^,-tj) 2 . The interpolated values resulting 
from the pairs (<*, tj), and (£, tj) 3 are used as left and right 
st ate* for the computation of one flux and (|,Jj ) 2 and 


Although the scheme just developed is based on the 
work of Levy, it is not to be confused with the algorithms 
of Ref. [5,6]. That work incorporates different interpolation 
methods, rotation angle selections, and integration schemes. 
Moreover, because of the current interpolation scheme, the 
cell edge formulation will not revert to a grid aligned 
scheme for arbitrary grids. The current formulation is used 
solely for simplicity and is probably not the optimal 
combination for a robust algorithm. These ’ingredients' of 
the recipe are selected so that maximum commonalty 
between the strategies is achieved. 

4.1 Cell-Center Orthogo nal Flux Implementation 

The imp lemen tation of the CCO strategy follows the 
same general outline as the previous CEO strategy. 
However, the rotation is performed at the cell center. First, 
the gradient at point (ij) is computed by the following 
stencil, 

bright = + 29i+l,J + 9i+l,;+l)» 

eft = -jfa.-U-l + 2 *~ 1.7 + *-1.7+1 ) » (4-23) 

<7up = ^(ft-l.y+l + 2q f iJ + 1 + ?i+l,;+l)» 

tfdown = J(«WJ-1 + 2 *.7-l + *+1.7-1 ) • 

The derivative in computational space is given by, 

4l - ( 424 > 

dq _ _ 

, dup 9down * 

Ol ? 

The gradient in physical space is given by Eq. (4.8). The 
metrics are taken as cell averages. The rotation angle is 
reputed using the same exponential blending function of 
Eqs. (4.9-4.11). However, a particular grid face must be 
se lected as the default direction ( 6 f ) which for this work, is 

the ( i+1/2 ) face. 


(£, i )) 4 are used for the other flux. 

The projection stage is computed by assigning the 
vectors of Eq. (4.4) as, 

n 1 =cos(0)i +sin(0)j , (4.21) 

n 2 = -sin( 0 )t~ + cos( 0 ); . 


Four interpolation points are required for the flux 
computation. The fluxes are computed with the cell center 
as one gtate and the appropriate interpolated point as the 
other grate The interpolation is computed using the same 
mapping methodology given in the previous section. In 
this case, however, the mapping is onto (£ 77 ) space where 
- 1£<*£1 and -l^rjil. 
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Before the fluxes can be projected, the algorithm must 
properly label each of the four fluxes. This sorting is 
required because the preferred direction may be oriented in 
any manner with respect to the computational coordinates 
(i.e. grid indices). So to project the "correct" fluxes onto 
the cell faces, the fluxes must be ordered. Flux sorting is 
accomplished by denoting the first flux counter-clockwise 

from cell normal ( i+1/2 ) as the n 1 direction (see Fig. 2). 
The al g o rithm selects the proper direction by the following 
criteria, 

n l • n M/2 ,j > 0 . n 1 x n i+U2J > 0 . (4.25) 

After the n 1 direction is defined, the other fluxes are 
incremented in the counter-clockwise direction. 

The last stage of the flux computation is the projection 
of the fluxes onto the grid faces. In the CEO method, the 
task is straightforward since the rotated fluxes are computed 
at each cell face and projected onto that face. However for 
the CCO method, the set of four fluxes must be projected 
around the entire cell. Such a projection is not unique so 
two methods that are implemented will be outlined. 

In the first method, an orthogonal pair of fluxes is 
assigned to each cell edge. With the notation 

F k =(F*n*)n*. (*=1,4), the fluxes at the cell edge are 

given as (see Fig. 2), 

^i+l/2,; = (^1 + £»)* "i+1/2,; - 

Kj+U2 = (^2 + ^1 ) * ".',,+1/2 • ( 4 - 26 ^ 

^,'-1/2 ,; = (^3 + F 2 ) * "i- 1/2,; • 

Kj-U2 ~ (^4 + * "i.;-l/2 * 

Recall that the fluxes are orthogonal so that the Jacobian is 
uni ty and the contravariant and covariant vectors are 
identical. 

The second method is a higher order projection 
(designated by HOP) where the flux across the grid face is 
not considered constant. In this case, the projected fluxes 
are weighted according to the local geometry. For example, 
consider the cell of Fig. 7. In computing the flux for face 
(i+l/2,j), the previous projection method (Eq. (4.26)) 
selects the flux pair F x and F 4 . The higher order 
projection divides the (i+1/2 j) face into two regions with 
the n 1 ray acting as the demarcation line. Then F 4 is 

projected onto the region below n} and F 2 is projected onto 
the upper region. The other component of the flux is 
obtained from the projection of F x over the region above 

the n 4 ray and F 3 below the ray. Another way to view 
this is that each flux acts on the cell faces or parts of cell 
faces in its half plane. The face fluxes are then given by, 


Fj+1/2 j = (fliFj + (l - q )Fj + 02 F 2 + (l - flj )F 4 ) • n M/2 j , 

Fi ji . U2 = (biF l +(l- bi)?, + b 2 F 2 +(l- «i.y+i/2 . 

(4.27) 

F;_,/2, ; - = ((l-ai)F, + <JiF 3 + (l- Oj )F 2 + OjF;) • «|_l/2.; 

Kj-m + 

where a x , h,, and b 2 are the weight coefficients. 
When computing the weight coefficients, the cell is taken 
to be a parallelogram because averaging the metrics over the 
cell, as is done for the interpolation, in effect reduces the 
actual cell to a parallelogram. Or in other words, the 
approximation in considering the cell to be a parallelogram 
in computing the weight coefficients is consistent with 
previous approximations. The projection is done over the 
actual cell to maintain conservation (in Eq. 4.27, the length 
scale is contained in the cell edge normal). The weight 
coefficients are computed by the following algorithm: first 
the coordinates of the dividing axes are given by. 


^l^cos^+^sin^), 
Vi = Vx cos(0) + T\ y sin(0) , 


= -|x sin(0) + $ y cos(0), 
772 = -Vx sin ( d ) + if; cos(0). 

Then, the switch functions are computed. 


max(&,T; t ) 


(k = 1 , 2 ), 


(4.28) 

(4.29) 


(4.30) 


Pu 


1 for|^|s|Tj t | 
0 for |6t| <M 


(* = 1,2), 


(4.31) 


S x = 
S 2 = 


1 for rjj 2 0 
-1 for rji < 0 

’ 1 for <j; 2 SO 
-1 for & < 0 


(4.32) 

(4.33) 


Finally, the weights are computed by, 

fli =1 — P 2 + yF2(l + r 2 ), (4.34) 

a 2 =i((l-F 1 )(l-5 1 ) + F 1 (l-5 1 ri)). 
b t =i(F 1 (l-S 2 )+(l-F l )(l-5 2 r2)). 

b 2 = P 2 + l 2 {l~P 2 ){l + r x ). 

Equations (4.28-4.34) may appear complex but are 
necessary to compute the weight coefficients without using 
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IF-THEN logic which inhibits vectorization and nearly 
doubles the amount of required computer time. 

4.4 Cell-Center Non-Orthogonal Flux Implementation 

The interpolation and logic required in the CCNO 
strategy is greatly simplified over CCO because everything 
is performed in computational space where the orientation is 
always the same. Also, since up winding is performed in 
each direction, a full range of rotation is achieved for 0 ° £ 8 
Z 90°. 

The gradient is computed in computational space using 
Eqs. (4.23-4.24). The rotation angle is given as, 

dq_ 

<4 - 35) 

where m is the magnitude of the gradient in 

computational space. Equation (4.35) is taken from Ref. 
[71. Also it should be noted that the values of DMIN for 
the blending function of Eq. (4.35) are quite different from 
the values of the exponential blending of Eq. (4.10) 
although the effect is nearly the same. Also, the values of 
6 given from Eq. (4.35) are reflected into the first quadrant. 

The interpolation for the CCNO strategy is much less 
complex than the previous two strategies because the 
topology in computational space is always consistent. In 
the CEO and CCO strategies, logic must be written to 
determine where to interpolate in physical space (Eqs. 4.12- 
4.20). For instance, depending on the skewness of a cell, a 
single interpolant value can lie in one of three possible 
quadrants in computational space. Conversely, in the 
CCNO method, a single interpolant value is restricted to 
only one quadrant. Therefore, the indices of the 
interpolation support stencil can be hardwired into the code. 


" 3 = (**,.»» + sin(0))/ 

+ (^uw cos ( 0 ) + Vw sin(0) V 

« 4= (v.- 1 cos ^)~^ sin ( e ))‘‘ 

Pairs of these vectors will be selected as the contravariant 
directions to be used in the projection for the cell faces. As 
before in the CCO strategy, two methods of projection are 
employed. The first assigns a flux pair to each face. The 
assignment follows the same schedule of Eq. (4.26) 
although now the full transformation must be computed. 
The second is the higher order projection where the rotated 
flux is considered to act on the regions of the faces in its 
half plane. The weight coefficients of the higher order 
formulation are much simpler to compute in CCNO 
because of the symmetry in the computational plane. First, 

let {ij} denote the pair of fluxes in the n 1 and n y directions 
which will be used for the projection of Eq. (4.4). The 
{zjjpair establishes the vector base for the projection. 
Further, let [ij}*n g denote this projection. Then the cell 
face fluxes are given by, 

F i+m .j = [^(^{1,2} + (1- W){3,4}) 

+a,{l,4}].n I+ i/ 2 ; . , (4.37) 

Ki*u 2 » [<h{W{ 2.3} + (1- 5W){1.4}) 
+a,{U}]«n i ,; +1/2 , 

Fi-mj = + (1 - SW){1,2}) 

+a,{2,3}]«n i _, /2i y , 
hj-m =[a2(5W{1.4} + (l-5W){2.3}) 

+a l {3.4}] • . 


The main difference in CCNO as compared to CCO is 
the flux projection. Recall that the fluxes are orthogonal in 
computational space but are non-orthogonal in physical 
space. Consequently, the projection of Eq. (4.4) does not 
reduce to simple trigonometric identities as in Eq. (4.22) for 
an arbitrary cell shape. However, this formulation reverts 
to a grid aligned scheme for arbitrary grids. First the 
following vectors are defined. 


where. 


SW = 


1 for <5 S r; 
0 for £ < f? 


a = min(£. 77 ) 


n,=i(l+a) , a 2 =j(l-a) . 


(4.38) 

(4.39) 

(4.40) 


cos (0)+ sin (6))i 

+ (^ cos ( 0 ) + ^ sin ( 0 ))^ 

i!= (v. cos(fl) ‘^ sin( 4 : 

+ (VU2 CO$ ( 0 ) " ZjH*, SiR ( 9 ty (4 ‘ 36) 


This completes the implementation of the CCNO 
method. 

4.5 Cell-Edge Non-Ortho gonal Flux Implementation 
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The CENO strategy is directly analogous to the 
previous CCNO method. The rotation angle calculation and 
interpolation are performed in computational space defined 



by the average metrics of Eq. (4.7). The rotation angle 
blending is given by Eq. (4.35) and the interpolating rays 
are given by Eq. (4.12) and Eq. (4.20). The projection onto 
the grid face is based on the vectors, 

n 1 = (f x cos(0)+ rf x sin(0))i + (|, cos(0) + rf y sin(0))> , 

(4.41) 

n 2 = (rf x cos(0) - ? x sin(0))r + (v y cos (0) - sin(0))) . 

The projection in the CENO method reduces to a simple 
trigonometric relation because the local metrics are 
constant. For example, for a 0+2/ 2) face the projection is, 

F; +1/2J = cos(0)- F 2 sin(0) . (4.42) 

The CENO formulation also suggests the possibility of 
selecting a preferred direction for the secondary flux as well 
as the primary flux (as opposed to setting a condition of 
orthogonality in physical or computational space). For 
instance, n 1 could be arranged to upwind across shock 
waves while n 2 is triggered to upwind across a shear wave. 
Such a scheme would require a full transformation of the 
fluxes and Eq. (4.42) would not apply. 

This completes the implementation of the four 
strategies. The next two subsections will briefly describe 
the boundary conditions and integration. 

4 6 Bound ary Conditions 

For a supersonic inviscid calculation, three types of 
boundary conditions are employed. The first is supersonic 
inflow where the flux is set to freestream on the inflow 
plane. The second is supersonic outflow where the 
variables are extrapolated to zeroth -order accuracy. The third 
is a nonporous wall condition where properties are reflected 
about the slip wall 

In the cell center strategy, the fluxes on the four sides 
of the cell are coupled to the rotation angle and hence to 
each other. Consequently, it is unclear how to define the 
rotation angle on the first row of cells next to the wall. So 
for this study, the fluxes on and near the wall are computed 
in a grid aligned fashion. More precisely, the nonporous 
wall fin** 1 * are computed based on the boundary condition. 
The flux on the faces normal and immediately adjacent to 
the wall are computed using a grid aligned formulation. All 
other fluxes are computed using the rotated flux functions of 
each particular strategy. This includes the outflow fluxes 
where a column of ghost cells is created with the flow 
quantities extrapolated to zeroth-order from the interior of 
the d omain This additional column of cells provides 
support for the interpolation stencil at the last column of 
interior cells. 


4.7 Integrat ion Method 

All four methods will be integrated in the same 
fashion. Thus, the only difference will be in the calculation 
of the cell residuals (RHS). The solution will be relaxed in 
time by the LU-SGS scheme of Yoon and Jameson 10 . 
Using the LU-SGS scheme with a rotated RHS is thought 
to be the best method to relax the solution to steady-state 
because of the following reasons: An explicit update 
scheme such as a multistep Lax-Wendroff type or Runge- 
Kutta requires multiple calls to the flux computation which 
should be avoided since the rotation requires extra work. If 
one were to seek an implicit solution, the temporal 
linearization of the fluxes would require some spatial 
approximation based on the local rotation angle since the 
fluxes are functions of interpolated quantities. This would 
expand the implicit stencil to a nonadiagonal block matrix 
for a two dimensional solution. With such a large banded 
matrix, an approximate factorization or Gauss-Seidel 
inversion becomes prohibitively expensive. Thus, the LU- 
SGS scheme is chosen. It is hoped that the diagonal 
Hnminanrp of the system, upon which the LU-SGS scheme 
depends, is sufficient to mask the incompatibility of the 
LHS to the RHS. 

5 0 Results 

The four rotation strategies are tested by computing the 
flow in a channel with a 15° ramp. The test case is taken 
from Ref. [5] where a more complete description is given. 
The 61x21 grid is shown in Fig. 8. The inflow Mach 
number is 2 and the flow is from left to right The test case 
is first computed to first-order accuracy using a grid aligned 
algorithm. Mach number contours are shown in Fig. 9. The 
range of the contours for this case and most subsequent 
cases are from 1.3 to 2.0 in increments of 0.05. The first- 
order grid aligned result shows a smeared initial shock wave 
and severely smeared reflected shock waves. 

The grid aligned algorithm is then extended to second- 
order accuracy by MUSCL extrapolation of the primitive 
variables. The MINMOD limiter is employed to damp 
oscillations around the shock wave. The resulting Mach 
contours are shown in Fig. 10. The second-order result 
shows the shock waves to be much more clearly defined. 

The next set of results are solutions of the four 
strategies where the rotation angle 0 is set to 45°. For this 
rhannftl flow problem, 45° nearly aligns the flux function 
across all the shock and expansion waves. Thus it presents 
a good initial test case for the rotated schemes without 
introducing any complications from a rotation angle 
rai niiatinn However, a priori knowledge of the rotation 
angles will not be known in general so the calculation of a 
preferred direction remains an important issue. 

The CE O strategy is presented first. All rotation 
angles are set to 45° on all faces (other than the particular 
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boundary faces discussed earlier). The Mach contours are 
shown in Fig. 11. Note for this case the increment in Fig. 

12 is 0.1 for plotting clarity, thus there are half as many 
contour lines as the other figures. The solution is highly 
oscillatory and shows signs of odd-even decoupling which 
appear as alternating shock waves and expansion shock 
waves. As currently implemented. Roe’s scheme admits 
expansion shock waves as solutions to the Riemann 
problem in violation of the entropy condition. The case is 
recomputed with a minimum eigenvalue condition. The 
new result is shown in Fig. 12 where the increment is 
returned to 0.05. The entropy correction significantly 
reduces the amount of odd-even decoupling (in comparing 
Fig. 12 to 13 keep in mind the number of contour lines). 
However, unacceptable oscillations still exist in the 
solution. 

The previous results demonstrate the shock capturing 
ability of the rotated flux that is the goal of the method. In 
spite of the oscillations, one can see that the shock waves 
are captured as sharply as the second-order grid aligned 
solution. However, the oscillations do exist and must be 
addressed 

The 45° case is now computed using CCO. The first 
solution uses the projection schedule of Eq.(4.26). The 
resulting Mach contours are shown in Fig. 13. The results 
show significantly improved shock capturing ability over 
the grid aligned results. Furthermore, the solution shows 
only a mild overshoot at the initial shock wave. The cell- 
center strategy captures the shock waves as well as the cell- 
edge technique but without the oscillations. This result is 
due to the averaging of the fluxes at the cell edge. Recall 
that each cell is rotated and the fluxes are projected to the 
cell edge. Therefore, each cell edge flux receives a 
contribution from each adjacent cell center. The averaging 
acts as a smoothing agent to the fluxes thereby reducing the 
amount of oscillations. This effect is a fortuitous by- 
product of the cell-center strategy. Also apparent from Fig. 

13 is an indication of an expansion shock wave. This will 
be discussed below. 

The next solution of Fig. 14 is computed using the 
CCO strategy coupled with the higher order projection 
described in Section 4.2. The result is identical to the 
previous case to plotting accuracy. 

The CCNO strategy is now used to compute the 
solution for a constant rotation angle of 45°. Recall that 
the rotation in this strategy is performed in computational 
space. Therefore, the rotation angle in physical space is 
somewhat different than the 45° in computational space 
depending on the cell shape and aspect ratio. The result for 
the standard projection schedule is shown in Fig. 15. The 
result is most interesting. The scheme captures the shock 
waves quite accurately. The first reflected shock wave is 
captured significantly better than the second-order grid 
aligned result of Fig 10. However, the solution also shows 


a captured expansion shock wave. Again, this is because an 
entropy condition is not included in the flux function. The 
presence of the expansion shock can be viewed as a positive 
result. This is an indication that the Riemann solver is 
producing results consistent with the one-dimensional 
operator independent of the grid. Not shown is the higher 
order projection result which does not affect the solution to 
plotting accuracy as in the previous cell-center strategy. 

The results of the CENO strategy are shown in Fig. 
16. This solution is computed with the entropy fix to 
reduce the oscillations generated from the cell-edge scheme. 
Recall, in this strategy the rotation is performed in 
computational space defined by the local metrics therefore, 
the fluxes across the cell edge are non-orthogonal. It is seen 
that the oscillations generated in this strategy are somewhat 
reduced from the CEO result (Fig. 12). However, the 
results are not as clean as either of the cell-center strategies. 

The convergence histories of the four strategies along 
with the first- and second-order grid aligned schemes are 
shown in Fig. 17. All schemes are run at the optimal CFL 
(the optimal CFL is determined by trial-and-error and is 
usually the maximum allowable time step but not always). 
The CFL numbers, and run time performance numbers for 
the 45° case are shown in Table 1. All solutions are run on 
a Cray-YMP and convergence is declared when the mass 
equation L2-Norm reaches 10'**. The "Time" column in 
Table 1 represents total run time including overhead. Since 
the run time of the first-order grid aligned computation is 
only 2. seconds, the time per grid point per iteration 
measure is somewhat misleading because the overhead is a 
significant percentage of the run time. 


Table 1 

Code Performance Parameters, 0= 45° 



CFL 

Iterations 

Time- 

seconds 

|isec/pt-iter 

(XI) GA 

oo 

248 

2.00 1 

6.72 

0(2) GA 

oo 

705 

5.14 

6.07 

CEO 

3. 

1276 

15.66 

10.22 

CCO 

4. 

715 

7.60 

8.85 

CCO-HOP 

8 

634 

7.19 

9.45 

CCNO 

8. 

1177 

11.41 

8.08 

CCNO-HOP 

8. 

1165 

12.70 

9.08 

CENO 

6. 

1304 

13.51 

8.63 


It is readily apparent that the rotation impedes 
convergence. This is caused by two effects. The first is the 
reduced dissipation and possibly increased dispersion from 
the rotation. The second is that the LHS does not account 
for the rotation on the RHS. Recall that all strategies 
employ the same LU-SGS LHS. It is also seen that the 
cell-center schemes converge faster than the cell-edge 
schemes. This is not surprising in view of the highly 
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oscillatory cell-edge results. The fastest converging rotated 
scheme is shown to be CCO-HOP. 

The next set of results is computed using the four 
strategies with the rotation angles based on the pressure 
gradient The angles are computed every time step and then 
frozen at a particular level of convergence or maximum 
iteration limit 

The CEO result is shown in Fig. 18. The solution 
based on pressure gradient does not show the signs of odd- 
even decoupling that were apparent in the 45° case. 
However, the solution is highly dispersive. The rotation 
angles are frozen after 300 iterations with a value of 
DMIN=Q.OOOOl. The frozen directions for the vertical faces 
are shown in Fig. 19 (every other line is plotted in the 
streamwise direction while every line is plotted in the 
crossflow direction). While viewing Fig. 19, keep in mind 
th^t an upwind flux is computed in both the preferred and 
normal direction. The horizontal face directions are similar 
to those presented in Fig. 19. The horizontal vectors next 
to the wall indicate that a grid aligned flux is computed. 
The directions are shown to basically align with the flow 
field shock waves. However, around the initial shock wave, 
the angles are shown to be oscillatory. The oscillations in 
the solution feed back to the rotation angle computation and 
create further oscillations in the solution. This effect can be 
reduced by filtering the rotation angle data as reported in 
Ref. [5]. 

The CCO results are presented next and are shown in 
Figs. 20 and 21. It is seen that the pressure gradient results 
are inferior to the 45° case but are still comparable to the 
second-order grid aligned solution. More precisely, the 
initial shock wave is captured better but the improvement 
Hitninichpt as the shock waves are reflected down the 
t-hamwl The rotation angles are frozen when the density 
residual reached a value of 5.x 10'^ and are shown in Fig. 
21. For this case DM IN=0. 0005. The directions are shown 
to align with the flow features and are well behaved. Recall 
that the angle computation algorithm is used in this 
case as in the previous CEO case. This is evidence that the 
ill-behaved directions of the CEO result are triggered by the 
cell-edge scheme. Again, no improvement is gained by the 
more complicated projection. 

The next set of results are for the CCNO method. The 
Mach contours are shown in Fig. 22. It is seen that the 
results are similar to the CCO strategy although slightly 
more dissipative. This is best seen in the first reflected 
shock wave. Also it is observed that the additional 
Higgipatinn , resulting from upwinding across the computed 
pressure gradient as opposed to setting the upwind direction 
to 45°, triggers the spreading of the expansion fan thereby 
eliminating the expansion shock wave. The rotation angles 
were frozen at the same convergence level as the previous 
CCO solution. A value of DMIN= 0.001 is used in 
conjunction with Eq. [4.35]. The rotation directions are 


mapped back into physical space and are shown in Fig. 23. 
The directions are shown to vary smoothly. Once again 
there is no improvement using the higher order projection. 

The last case computed is the CENO strategy with the 
rotation based on pressure gradient The result is shown in 
Fig. 24. The solution is very similar to the CEO result of 
Fig. 18 although the non-orthogonal projection is slightly 
less dispersive. The rotation angles were frozen after 300 
iterations and are shown in Fig. 25 (DMIN= 0.001). 

The convergence histories of the rotated solutions based 
on pressure gradient are plotted in Fig. 26. The CFL 
numbers and run time data for the pressure gradient cases are 
shown in Table 2. The following observations are made. 
First, the cell-edge schemes do not converge until the 
rotation angles are frozen at 300 iterations. After the 
freezing, the convergence rate is less than the cell-center 
schemes. The cell-center schemes converge faster than the 
second-order grid aligned case and, interestingly enough, 
taifp. less CPU time to compute a comparable solution. 
From experience, the cell-center strategies will converge 
about 3 to 3.5 orders of magnitude before a limit cycle is 
reached. Freezing the angles allows complete convergence. 


Table 2 

Code Performance Parameters, 9 based on Vp 



CFL 

Iterations 

Time- 

seconds 

psec/pt-iter 

0(1) GA 

oo 

248 

2.00 

6.72 

0(2) GA 

oo 

705 

5.14 1 

6.07 

CEO 1 

3. 

789 

10.42 

11.01 

CCO 

20. 

376 

4.26 ! 

9.44 

CCO-HOP 

OO 

360 

4.33 

10.02 

CCNO 

20. 

389 

4.04 

8.65 

CCNO-HOP 

OO 

399 

4.58 

9.57 

CENO 

7. 

790 ”1 

8.44 

8.90 


Throughout this study, the cell-edge schemes are shown 
to consistently create significant oscillations in the solution 
of this channel problem. There have been a variety of 
approaches to reduce these oscillations. Levy computes the 
45° rotation case by only rotating the horizontal faces and 
computing a grid aligned flux on the vertical faces. This 
eliminates the oscillations but the improvement in shock 
capturing ability is reduced. In the pressure gradient case, 
averaging of the rotation angles throughout the domain is 
employed to reduce the feedback between the flow field 
oscillations and the rotation angle computation. Another 
technique is the simplified interpolation scheme of Ref. [6]. 
The cim piifirel scheme reduces the dependence of the flux on 
the rotation angle since small changes in the angle do not 
change the interpolated dependent states. The cell-center 
strategies, on the other hand, allow greater flexibility in the 
angle selection and the interpolation. This is because the 
flux averaging, required in order to define a unique flux at 
the face, acts as an inherent smoothing agent in the scheme. 
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6 0 Conclusion 

Four rotated upwind strategies have been tested on a 
simple Mach 2 channel flow problem* It is concluded that 
the cell-center strategies offer better promise in extending 
the study to three dimensions than the cell-edge strategies 
because of the reduction in work, the inherent smoothing in 
the flux averaging, and acceptable convergence rates. The 
CCO strategy produces the best results of the four 
strategies. The results are comparable to the second-order 
grid aligned scheme. The CCNO strategy also produces 
good results, is simpler to incorporate than CCO, and has 
the added feature of reverting to a grid aligned formulation. 
Finall y, it is shown that for this simple test case, the 
higher order projection does not improve the cell-center 
results. 
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Figure 1.) Cell-Edge Orthogonal Flux Strategy. 
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Figure 9.) Mach Contours for First-Order Grid 
Aligned, Increment = 0.05. 



Figure 10.) Mach Contours for Second-Order Grid 
Aligned, Increment = 0.05. 



Figure 11.) Mach Contours for CEO, #=45°, 
Increment = 0.1. 



Figure 12.) Mach Contours for CEO, #=45°, Entropy Fix, 
Increment = 0.05. 





Figure 13.) Mach Contours for CCO, &=45°, 
Increment = 0.05. 



Figure 14.) Mach Contours for CCO, Higher Order 
Projection, &=45°, Increment = 0.05. 



Figure 15.) Mach Contours for CCNO, 0=45°, 
Increment = 0.05. 



Figure 16.) Mach Contours for CENO, <M5°, Entropy 
Fix, Increment = 0.05. 
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Figure 17.) Convergence Histories for 9 = 45 ° Solutions. 



Figure 18.) Mach Contours for CEO, 9 based on Vp, 
Increment = 0.05. 



Figure 19.) Frozen Rotation Directions on the 
vertical faces for CEO, 9 based on Vp. 



Figure 20.) Mach Contours for CCO, 9 based on Vp, 
Increment = 0.05. 
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Figure 21.) Frozen Rotation Directions for CCO, Figure 24.) Mach Contours for CENO, 9 based on Vp, 

9 based on Vp. Increment = 0.05. 



Figure 22.) Mach Contours for CCNO, 9 based on Vp, Figure 25.) Frozen Rotation Directions on the vertical 

Increment = 0.05. faces for CENO, 9 based on Vp. 



Figure 23.) Frozen Rotation Directions for CCNO, 
Abased on Vp. 



Figure 26.) Convergence Histories for 9 based on Vp 
solutions. 
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Abs tra ct 

Rotated upwind algorithms are presented for the 
numerical solution of the Euler and Navier-Stokes 
Equations in two and three dimensions. The finite-volume 
algorithms are designed with the notion of aligning Roe’s 
approximate Reimann solver in a computed preferred 
direction. The algorithms are applied to shock wave 
reflection and shock-wave/boundary-layer interaction 
flowfields typically found in supersonic and hypersonic inlet 
configurations. Calculation of an inviscid Mach 2 channel 
flow problem shows that the rotated algorithm produces 
more accurate results than a traditional grid aligned 
algorithm to both first- and second-order accuracy. 
Moreover, the improvements to second-order accuracy are as 
great as those to first-order accuracy. Viscous solution of a 
turbulent boundary layer / shock wave impingement show 
that the rotated scheme improves the shock wave capturing 
in the inviscid portion of the flowfield to both first- and 
second-order accuracy. The improvements in the shock 
wave capturing to first-order accuracy result in improved 
wall pressure and skin friction distributions. However to 
second-order accuracy, the wall predictions are not improved. 
The calculation of an inviscid three-dimensional shock wave 
surface shows the rotated algorithm to be more accurate than 
the grid aligned algorithm to both first- and second- order 
accuracy. The accuracy improvements in three dimensions 
are not as great as those in two dimensions. Computation 
of a viscous flowfield in the corner of two intersecting 
wedges shows that the inviscid portions of the flowfield are 
qualitatively improved with the rotated algorithm to both 
first- and second-order accuracy. However, surface pressure 
predictions are only marginally improved with the rotated 
algorithm. 

1.0 Introduction 

Much of the recent algorithm development has been 
focused on creating artificial dissipation models that include 
multi -dimensional effects. It is hoped that the great success 
of the upwind dissipation models can be extended by more 
accurately modeling wave propagation in multi-dimensional 
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space. A variety of promising and novel approaches are 
currently being pursued. Since fhe number of methods has 
become too numerous to list in this format, the reader is 
referred to a comprehensive overview of the multi- 
dimensional methods presented by van Leer 1 . 

One approach to including multi-dimensional effects in 
the flux function is the use of a rotated upwind method. 
The idea, originating with formulations of the potential 
equations 2 ’ 3 and later developed by Davis 4 for the Euler 
equations, is to rotate the integration stencil to an 
orientation where the differencing is applied in a direction 
that more closely models the correct path of information 
propagation. For the Euler and Navier-Stokes equations, 
rotation results in applying the upwind formulation across 
captured discontinuities rather than along grid lines. The 
orientation of the upwind operator with respect to flow 
discontinuities significantly impacts the ability of the 
upwind model to correctly interpret local wave propagation. 
For example, Roe's^ scheme is known to model a single 
shock wave as a sum of acoustic and shear waves if the 
operator is applied oblique to the shock wave (an excellent 
discussion of this effect is presented in Ref. [6]). The 
additional non-physical shear wave generates excess 
dissipation and smears the captured shock wave. However, 
the correct wave propagation is predicted when the operator 
is applied across the shock wave. Thus, we are motivated 
to calculate an upwind flux across flow discontinuities 
independent of the orientation of the discontinuity with 
respect to the grid. 

Many rotated upwind schemes for solution of the Euler 
equations have been developed^’^lO^ all of which show 
increased shock resolution as compared to traditional grid 
aligned schemes for simple inviscid test problems to first- 
order accuracy. The improvements to second-order accuracy 
are shown to be modest but encouraging. Furthermore, 
none of the previous rotated upwind schemes have been 
extended to three dimensions. In fact, very few of the newly 
developed multi-dimensional upwind methods have been 
tested in three dimensions. Some exceptions are the works 
of Deconinck et. al. 11 for three-dimensional scalar 
advection, and Rumsey 6 who solves the Euler equations in 
three dimensions. Rumsey's results show improvement of 
three-dimensional shock and shear capturing to first-order 
accuracy. To second-order accuracy, only improvement to 
shear capturing is shown. Thus, it has not been 
demonstrated, nor is it necessarily anticipated, that a rotated 
upwind scheme that aligns the Roe operator in a single 
preferred direction in three-dimensional space will improve 



three-dimensional shock wave capturing. It is the objective 
of this study to develop a robust, rotated upwind algorithm 
for the solution of both the Euler and Navier-Stokes 
equations in two and three dimensions. The necessary 
characteristics of this scheme are that it is fully 
conservative, second-order accurate in space, and sufficiendy 
robust to compute a variety of supersonic through 
hypersonic flowfields. 

2 0 The Navier-Stokes Equations 

For this study, the governing equations are the set of 
conservation laws governing three-dimensional viscous fluid 
flow known collectively as the Navier-Stokes equations. 
The equation set in integral form and in index notation is 
given as the following: 

— lpdV+ J(pu,n 1 )dS = 0, (2.1) 

Volume Surface 

— J pujdV + J (pup j + Ty )n;dS = 0 , (2.2) 

^ Volwru Surfazt 

■j jpE,dV+ j(E l u i +T ij Uj+q i )n i dS = 0, (2.3) 

Volume Surface 

where p is the density, is the velocity, E t is the total 
energy, is the outward normal to the surface bounding 
the control volume, is the heat flux, dV is the 
differential element of the control volume, and dS is the 
differential element of the bounding surface. The term T tj 
is the stress tensor which contains both the pressure and 
shear stress terms. The shear terms are modeled assuming a 
linear stress-strain relation, the heat flux is modeled through 
Fourier’s law, and the dynamic viscosity is calculated using 
Sutherland’s Law. In addition, the perfect gas equation of 
state is employed. 

The governing equations are solved in a finite-volume 
formulation with the conserved variables taken as averages 
over the cell volume. The surface integrals of Eq. (2.1 - 
2.3) become summations over the cell sides. Expanding the 
vector momentum equation into its scalar components, the 
equations for a fixed grid become, 

i^ = I — Y F»S, (2.5) 

dt Volume 

where U is the vector of dependent, conserved variables and 
F contains the fluxes. The cell face normal as expressed in 
terms of the grid metrics is denoted by S . The magnitude 
of S is equal to the cell face area. The Euler equation 
subset is obtained by neglecting shear terms. 


I n Two-Dimensional Rotated Upwind Algorithm 

Traditional upwind schemes solve the Riemann 
solution or split the flux Jacobian along grid lines. Or, in 
other words, the upwinding is performed in the grid 
contravariant directions regardless of the orientation of the 
flow features with respect to the grid. An alternative to the 
grid aligned method is the rotated upwind scheme which 
aligns the upwinding stencil in a direction that is more 
likely to model the flow physics. This section describes 
such a rotated upwind algorithm for the solution of the 
governing equations in two dimensions. This is followed 
by Section 4.0 which presents the extension of the 
algorithm to three dimensions. 

The two-dimensional algorithm is an extension of the 
first-order Cell-Centered Non -Orthogonal flux (CCNO) 
strategy presented by Kontinos and McRae The 
approach is to compute a rotation angle in computational 
space based on the gradient of a flow field scalar such as 
pressure or Mach number. The difference stencil is then 
rotated into the preferred orientation as is shown in Fig. la. 
The rotated axes define four rotated contravarient directions 
which are functions of the cell face metrics and the rotation 
angle. The primitive variables are then interpolated onto 
the new stencil through linear interpolation between 
adjoining cell centers. The four fluxes are computed in the 
rotated contravarient directions using Roe’s scheme with the 
cell center as one state and the appropriate interpolation 
point as the second state value. The next step is to project 
the fluxes back onto the original grid faces. However, the 
orthogonal axes in computational space become four non- 
orthogonal directions in physical space as shown in Fig lb. 
Therefore, in order to project the fluxes in a conservative 
manner, each cell face projection must be computed using a 
full metric tranformation. This results in four metric 
tranformations being computed at every cell for every 
iteration. The rotation angles are computed at every cell 
center every iteration for 2-3 orders of magnitude 
convergence and are subsequently frozen. Complete details 
of this algorithm can be found in Ref. [12, 13]. 

The CCNO algorithm is extended to second-order 
accuracy using the MUSCL philosophy wherin the values 
of the state variables are extrapolated to the interface. In a 
standard grid aligned formulation, the extrapolation is 
performed along grid lines. For example, the right state of 
the (/+1/2, j) face uses point (i+2, j) as support for the 
extrapolation. Similarly, the right state of the (i,j+ 1/2) 
face uses point (t, j+2) as support. In the CCNO 
form niati nn, the same procedure is implemented, however, 
the second-order support point is a function of the rotation 
angle. For example, consider the first quadrant of the 
rotated stencil in computational space as shown in Fig. 2. 
The first-order line of interpolation is represented as the 
dotted line in the figure. Also shown is the next ring of 
grid points around point (i, j) which provides support for 
the second-order interpolation. Based on the rotation angle. 
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one could follow the same liner interpolation procedures as 
the First-order stencil. However, the second-order support 
ring consists of six points thus requiring logic to determine 
the location of the interpolation. Such a procedure 
significantly increases the amount of work required by the 
algorithm. An alternative procedure is implemented in this 
study. The second-order support point is taken to be the 
closest cell center value to the intersection of the rotated 
axis with the second-order support ring. Thus, the second- 
order interpolation is not a continuous function of the 
rotation angle, but is a step function. Of course, the 
introduction of a discontinuous interpolating function 
dictates that the rotation angles must be frozen in order to 
reach convergence. If the angles are not frozen, it is very 
probable that the second order interpolation will cycle 
between two discontinuous values. But since it is seen in 
Ref. [12] that the rotation angles are frozen in order to 
converge the first-order algorithm, the discontinuous second- 
order interpolating function does not introduce any further 
limitations. 

The Final component of the flux computation of the 
CCNO algorithm are rotated boundary conditions. Since 
the rotation angles are defined at the cell centers, the flux 
computation of each of the four faces surrounding a 
particular cell center are essentially coupled. Consequently, 
reflection and extrapolation boundary procedures commonly 
employed in grid aligned schemes must be modified to 
account for the coupling of the boundary faces with 
adjoining interior faces. The details of the boundary 
condition procedures are omitted here for brevity but are 
presented in Ref. [13]. Finally, the viscous fluxes are 
computed in the standard grid aligned fashion and the 
solution is relaxed in time using the diagonal form of the 
LU-SGS scheme of Yoon and Jameson 14 . 

4 0 Three-Dimensional Rotated I Tnwind Algorithm 

This section will oudine the development of the three- 
dimensional cell-centered non-orthogonal flux rotated 
upwind scheme. Recall, that the overall strategy of the 
CCNO scheme is to compute a preferred upwind direction in 
computational space at each cell center. Then a coordinate 
rotation is performed such that a computational axis is 
aligned in the preferred direction. Primitive variables are 
then interpolated onto the new coordinate system. The 
inviscid fluxes are computed and transformed back onto the 
grid faces in a conservative manner. Sections 4. 1-4.4 detail 
the steps of the algorithm. Section 4.1 presents the 
coordinate rotation that aligns a coordinate axis with any 
arbitrary direction. The rotation matrix deFined in Section 
4.1 is then used in Section 4.2 to define the contravariant 
directions of the rotated fluxes. Section 4.3 presents the 
interpolation procedure and Section 4.4 describes the 
projection of the rotated fluxes onto the original grid faces. 


4 1 Coordinate Rotation 

The primary task of aligning a computational axis with 
an arbitrary direction can be performed in a variety of ways. 
For instance, given some direction, three Eulerian angles 
can be found which can be used to set up a rotation matrix. 
Application of this rotation matrix will cause one of the 
original coordinate axes to align with the preferred direction. 
However, this method of rotation will align the same 
original axis in the preferred direction for every given 
orientation. In other words, after rotation the original 
coordinate axis can end up in any of the eight quadrants of 
the original system. So upon calling the interpolation 
routine, which is required to compute the fluxes, logic must 
be included to determine the orientation of the new system 
with respect to the old. Therefore, the secondary objective 
is to minimize the number of quadrants in which a 
particular axis can lie. In a practical sense, the attempt is to 
define a rotation such that the indices of the interpolating 
stencil can be hardwired into the code, thereby reducing the 
amount of logic and computer work. The objectives are 
achieved by defining two possible rotational 
transformations. Application of these two transformations 
is g uar anteed to align some computational axis with any 
given preferred direction. The following discussion 
describes the two rotational transformations and their 
implementation into the algorithm. This development is 
aided by the work of Mayer 15 who presents an excellent 
discussion of rotation matrices. 

The first rotational transformation, Type I, is 
characterized by a rotation through B about the C axis 
denoted by Z(9) then a rotation about the <£;' axis through 0 
denoted by 2(0). The rotation is restricted such that 

0 < 6, 0 < 90° . The transformation is then expressed by, 

(£' V' £') = (« *1 C)Z( 0)2(0). (4.1) 

Figure 3 depicts the Type I rotation. The comers of the 
cube are cell centers in computational space. The center of 
the cube is point For clarity, not all of the cell 

centers of the interpolating stencil are drawn. Notice in 
which quadrant of the original system the new axes lie. The 
+ 77 ' axis lies in the (-,+,+) quadrant, -rj':(+ +£':( + . -■,+)> 
So given a preferred direction in any of these 
four quadrants, the Type I rotation will align one of the 
computational axes given a proper definition of the rotation 
angles. Also note that the +£' axis lies on the border of the 
(+,+,+) quadrant and 

As can be expected, the Type II rotation will allow 
alignment of an axis given a preferred direction in the 
remaining four quadrants. The Type II rotation is given by 
a rotation through 6 about the f axis then a rotation about 
the r\' axis through -0 denoted by H(-0). The negative sign 
on 0 is a consequence of the right-handedness of the system. 
This transformation is expressed by, 
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U' v' C)=(Z n f)z(0)H H) . (4.2) 

Figure 4 depicts the Type II rotation. Now it is seen that 
the +£' axis lies in die (+,+,+) quadrant, -? 
+£•;(+,-,-), The +r?' axis lies on the border of the 

(-,+,+) quadrant and Comparing Type I and Type 

II rotations it is seen that the ends of <jj' and tj' always lie in 
the same quadrant. Thus, for the and Tj axes the 
objective of minimizing the number of possible 
orientations has been achieved. It is only C which requires 
some logic to determine indices of the interpolating stencil. 
However, each end of the £' axis can only he in one of two 
quadrants. Superposition of the two transformadons result 
in the following rotadon matrix. 




1 for pr > 0 

» 

[0 for pr < 0 


(4.7) 


S^aS^+pS^. (4.8) 

Then from geometric considerauons the angles are given by, 


cos( 0 ) = - 


1 


2 2 
P +4 


(4.9) 


R(a,9, </>) = 


cos (6){a + P cos{(t>)) 
-sm(9)[P + a cos( 0 )) 


cos (<t>) = 


(l - S * )|r| + 5* Vp 2 + 9 2 

Vp +<7 2 + r2 


(4.10) 


|_sin( 0 )(asin( 0 )-/Jcos( 0 )) 
sin(9)(a + P cos(<p)) Psin{ip) 

cos{9)(P + acos( 0 )) asin( 0 ) , 

-sin(^)(acos( 0 ) + /Jsin( 0 )) cos (<t>) _ 


(4.3) 


where is a = 1 for a Type I rotadon, a = 0 for a Type II 
rotation, and f5 = 1 -a 

To completely define the rotadon, the angles and the 
switch between Type I and Type II rotadon must be 
computed. First let the preferred direcuon in computadonal 

space be given by =^p| + < 7 T + r C) where the A 

denotes some unit normal in computadonal space. The 
switch between types of rotations is given as. 


The complete range of rotation is 0 < 9,<t> < 90° so that the 
cosine of the angles is sufficient 

4 2 Contra variant Directions 

The rotation matrix developed in the previous section is 
used to define the contravariant directions of the new 
coordinate system. Recall that the advantage of defining the 
rotadon in computational space is that the scheme will 
automatically revert to a grid aligned formulation. 
Moreover, the scheme reverts to a grid aligned mode in the 
limit of 90° rotation. Conceptually, the coordinate 
transformation rotates each coordinate axis between 
adjoining cell face normals which are the contravariant 
directions of the grid. For example, let J be the metric 
Jacobian matrix given by. 


1 for pq <0 




0 for pq'Z 0 


(4.4) 


/ = 


Vx Vy ^ 1 : 

c, c, c. 


(4.11) 


Before computing the rotation angles, the orientation of the 
preferred direction is revealed by computing the following 
switch functions. 


In a finite-difference formulation, the metrics are defined at 
the grid point so that the metrics of a new coordinate 
system, denoted by J' are unambiguously given by. 


S 


pi 


f 1 for pq> 0 
[0 for pq <0 


f 1 for qr> 0 
[0 for qr < 0 


(4.5) 


(4.6) 


J' = R(a,9,(l>)J . (4.12) 


In a finite-volume formulation, the metrics are defined at the 
cell face. Therefore, under a cell centered coordinate 
rotation, the new metric Jacobian is not clearly defined. 
However, the rotation matrix does prescribe the proper 
transition between the cell face normis under a rotation in 
computational space. Consider the +£ rotated contravariant 
direction as an example. We know that under the given 
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rotation, the +£' axis always lies in the first quadrant of 
computational space. Furthermore from observing the 
limits of 0 ° or 90° rotation of both 8 and 0 , it is apparent 
that the contravariant direction is some function of the 
metrics of faces (i+1/2, j, k ), (i, 7+1/2, k), and (i,7* £+1/2). 
Therefore the transformation of the +£' axis is given by, 



^*•♦1 fXj* 

n.j± 

^ z i+in.i4 

CD 

c< 

II 

Vx 


Vz 



Cy 

'LiMin 

1 

C 

vjl 


. (4.13) 


It is important to note that this transformation yields the 
direction of the +;' axis only, consequently, only the top 
row of elements of J ' is needed giving, 

n x = j[ x i + j[ 2 j + j\-$k , (4.14) 

where 7 ^* are the elements of /' and n 1 is the contravariant 
direction of the +g coordinate axis. 

The other five contravariant directions are defined by 
developing some simplifying notation. First, let J{ci t b t c) be 
defined as follows, 





a 

1 +a.M 

J(a,b,c) = 

Vx. 

77 y 



J 

1 


c, u 

1 


Here, a , Z>, and c assume values of ±1/2 to select the sides 
of the computational cube at point (i,j, k). Furthermore, 
let the function ROW(m, A) give the vector on the mth 
row of matrix A . Then the n x direction developed 
previously is denoted as 

n 1 = ROW(l, R(a,9,<t>)J(ll2,U2,U2)) (4.16) 

The remaining contravariant directions are developed in a 
similar manner and are given as, 

n 2 = ROW(2, R(a,6,4>) /(-t 72 , 1 / 2 , 1 / 2 )) , (4.17) 

M 3 = ROW(l, R(a,9,<p) J{-\ / 2 , - 1 / 2 , - 1 / 2 )) . (4.18) 
n A = ROW(2, R(a,6,<j>) /( 1 / 2 ,- 1 / 2 ,- 1 / 2 )) . (4.19) 


n 5 = ROW(3, R(a,9,<t>)J{a-u2,-\i2,ii2)) , (4.20) 


n 6 = ROW(3, R{a, 0 , 0 ) 7 (i / 2 - a,i / 2,-1 / 2 )) . (4.21) 


where n 2 , n 3 , n 4 , n 5 , and n 6 are the contravariant 

directions of the -£\ +tj‘, -Jj\ + C> ~C axes, respectively. 


4 3 Interpolation 

Once the rotation of Section 4.1 has been defined, the 
interpolation of the primitive variables can be performed. 
In fact, one of the advantages of rotating in computational 
space is that there is symmetry of the interpolation stencil 
about the cell center. The first-order interpolation is 
performed by computing coordinates of the interpolating ray 
and then using a trilinear interpolating function over a 
quadrant of computational space. For example, consider the 
( 4 .'+'+) quadrant of Fig. 5 where the n * ray is shown under 
a Type II rotation. Let the origin of (£, tj, Q be at point 
(j,;, k). Then over the region of the cube any value q can 
be given by the trilinear function, 

q = (1- C)[fr i+ i.M + Vq iJ+ i.t -(!-?)(!- 

+ qv(qi+\J+l,k ~ ?i+l ,j.k ~ 4iJ+ l,t )] (4.22) 

+ c[^i+l,/.t+l + 7 7 < 7i,/+l.t+l — (l - ?)(^ — v)qi,j.k + 1 

+ ^ T /( < 7.+l,y+l.t+l _ q i+\.j ,k+\ - 9i,/+l,*+l)]‘ 

All that needs to be defined are the coordinates Ojj.rj.f) in the 
range 0 £ tj,£ 5 1. Of course it is desirable to define the 
coordinates such that for no rotation, a grid aligned scheme 
is recovered. One possible option is to set (^. p, 0 such that 

•^ 2 + tj 2 + C 2 =1. This is shown as the shell segment 

drawn in Fig. 5. This is the easiest option to implement 
because all that is required is a normalization of the 
coordinates defining the interpolating ray. The option used 
in this study is to clip (£,77 . 0 where the interpolating ray 
intersects the outer face of the interpolating stencil. Then 
the trilinear function of Eq. (4.22) in effect reduces to 
bilinear interpolation over the outer square region. 
Moreover, if the interpolating ray happens to intersect the 
line segment adjoining two grid points, the interpolating 
function reduces to linear interpolation. This second option 
is selected so that under/overshoots are minimized through a 
low order interpolating function. Similar expressions to 
Eq. (4.22) are written for the other quadrants of the 
interpolating stencil. 

The second-order interpolation of the three-dimensional 
algorithm is similar to the two-dimensional procedure of 
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Section 3. The second-order interpolant value is taken as 
the nearest cell center value to the interpolating ray. 
Consider the first quadrant of the second-order stencil as 
shown in Fig. 6. In the first quadrant, there are 19 points 
in the stencil which are displayed in Fig. 6 as black and 
gray circles (excluding (i, j, *)). The discontinuous 
interpolation function must be able to select the nearest cell 
center value to the interpolating ray. Such a function is 
costly because there are 19 possible interpolant values. 
Consequently, the stencil is reduced to be composed of the 7 
comer points which are shaded black in Fig. 6. The indices 
of the nearest cell centers are calculated based on the 
coordinates of the first-order interpolating ray. Details of 
this procedure are contained in Ref [13]. 

4 4 Projection 

The inviscid fluxes are computed using Roe's scheme 
with the cell center point as one state, the appropriate 
interpolated point for the second state, and the rotated 
metrics defined in Section 4.2. At this point in the 
computation, there are six fluxes computed in the rotated 
contravariant directions. In order to maintain conservation, 
the fluxes must be projected back onto the original grid 
faces. Recall that the fluxes are orthogonal in 
computational space; however, they are non-orthogonal in 
physical space. Therefore, in order to project the fluxes, a 
new transformation must be computed for each face based 
on the contravariant directions as is done in the two- 
dimensional algorithm. The three-dimensional 
transformation can be found in Thompson et. al.^ . Once 
again for emphasis, the full transformation must be 
computed at each grid point for each face in order to 
maintain conservation. This results in six metric 
transformations being computed at each grid point every 
iteration. Unfortunately, this requirement causes significant 
increases in the algorithm's work load. 

The projection schedule for each of the six faces of the 
cell is now given. Three fluxes are required per face and are 
selected such that in the limit of 90° rotation the 
contravariant directions are properly oriented in a grid 
aligned fashion. For example, the F f • S M f2j 9 k ( w ^ ere 

F i is the inviscid flux and S,+i/ 2 ,),* * s norma ^ 

of face (i+i/2, j, k)) must be a function of the contravariant 
directions that have a (i+7/2, j, k ) metric dependency. From 
Eq. (4.16) and Eqs. (4.19 - 4.21) it is seen that n\ n 4 , n 5 , 
and n 6 all have a (i+7/2, j, k) metric dependency. Letting 
H represent the transformation, the projection dependency 
for face (i+7/2 J, k) is given by, 

VW'f'W!' <423) 

where = F j •n* and ct ^nd ^ sre the switching 
coefficients between Type I and Type II rotation. The 
projection schedule for the remaining faces are given as. 


•S i -u2j.k = H{F 2 'F 3 JF 5 .aF 6 ) , 

(4.24) 

Ft • j+\n.k — ) • 

(4.25) 


(4.26) 

F,'S iJMm = H(F\F 2 ,F 5 ) , 

(4.27) 

F/-S, M _ 1/a = //(F 3 ,F 4 ,F 6 ). 

(4.28) 


The two numerical flux values from each adjoining cell 
center rotation around a particular face are averaged in order 
to uniquely define the cell face flux. Also, rotated boundary 
conditions are implemented in three dimensions and are 
described in Ref. [13]. 

5.0 Results 

5.1 Inviscid Mach 2. C hannel Flow 

The two-dimensional CCNO algorithm is tested by 
computing a Mach 2 flowfield in a channel with a 15° ramp 
which can be found in Ref. [7], This geometry has been 
used quite extensively as a test case for multi-dimensional 
upwind algorithms. The deceptively simple geometry 
belies the intricate physics of the flow. In fact, fine grid 
resolution studies of this case by Parpia and Parikh 1 show 
that the initial shock wave does not reflect in a regular 
manner off of the top wall. Instead, a Mach stem develops 
and a slip surface is generated at a triple point. These flow 
characteristics are not captured with the 61x21 grid (Fig. 7) 
used in this study. The case is computed using the standard 
grid aligned upwind formulation and the resulting Mach 
contours are shown in Fig. 8. The contours range from 1.3 
to 2.0 in increments of 0.05. The case is computed using 
the CCNO algorithm to first-order accuracy. The preferred 
upwind direction is selected to be the direction of the 
pressure gradient. The resulting Mach contours are shown 
in Fig. 9. The rotation angles are computed for every cell 
every iteration until the L 2 -Norm of the mass equation is 
reduced by approximately three orders of magnitude 
whereupon the rotation angles are frozen. The freezing 
occurs in 260 iterations and the resulting orientations are 
shown in Fig. 10. The solution converges 12 orders of 
magnitude in 453 iterations at a CFL of 5 consuming 2.9 
seconds of Cray Y-MP/C90 time. 

The first-order CCNO method is shown to capture the 
shock waves with less smearing than the first-order grid 
aligned result of Fig. 8. The solution is relatively 
oscillation free except for a leading overshoot ahead of the 
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ramp shock wave. A pressure survey through the middle of 
the domain is shown in Fig. 11. The pressure comparison 
clearly shows that the rotated scheme captures the 
discontinuities and associated pressure rises more sharply 
than the grid aligned procedure to first-order accuracy. In 
order to quantify the error, a fine grid solution of 385x129 
points is computed using the second-order grid aligned 
scheme. This solution is considered to be the “exact” 
solution of an error norm calculation. Mach number 
contours of the 385x129 grid solutions are shown in Fig. 
12. Based on the fine grid solution, the density L 2 -error- 
norm of the first-order grid aligned solution is calculated to 
be 4.561xl0' 3 while the first-order CCNO error is 
2.521xl0' 3 . Analysis of the error distribution reveals that 
the major source of error is the smearing of the captured 
shock waves. 

The increase in accuracy of a first-order rotated scheme 
has been shown by previous researchers e.g. Levy et. al. 7 , 
Dadone and Grossman 8 , and Kontinos and McRae 9 - 12 . 
However, the improvements of second-order accurate rotated 
schemes over second-order accurate grid aligned schemes 
have been shown to be modest at best. Thus, the CCNO 
algorithm is now compared to grid aligned upwinding to 
second-order accuracy on the 61x21 grid. The Mach 
contours of the second-order grid aligned and CCNO 
solutions are shown in Fig. 13 and 14, respectively. The 
range and increment of the contours are the same as the 
previous figures for direct comparison. The frozen rotation 
angles, based on the pressure gradient, are similar to those 
in the first-order solution as shown in Fig 10. From 
comparison to the second-order grid aligned result, it is seen 
that the rotated solution shows improved shock wave 
capturing. The improvement in the solution is also 
indicated in the pressure survey along j = 10 shown in Fig. 
15. The density error norm of the grid aligned and CCNO 
solution are calculated to be 2.640x10’ 3 and 1.345x10 3 , 
respectively. Two very important observations can be made 
from the error calculation. First, the first-order CCNO 
solution is more accurate than the second-order grid aligned 
result. Secondly, and most importantly, the ratio of the 
error norms show that the improvement of the CCNO 
solution over the grid aligned solution is as great to second- 
order accuracy than to first-order accuracy. The CCNO 
solution with CFL=5 converges 12 orders of magnitude in 
897 iterations in 7.33 seconds on the Cray Y-MP/C90. 

The density error norms are plotted versus the number 
of points in the x direction on a log-log scale in Fig. 16. 
Additional data points, generated from 97x33 and 121x41 
grids, are also included in Fig. 16. The figure shows that 
the rotated results are more accurate than the grid aligned 
results to both first- and second-order accuracy. 
Furthermore, the improvements to second-order accuracy are 
as great as to first-order accuracy. It is also seen that the 
accuracy gains are maintained as the grid is refined. Also 
shown in Fig. 16 are the slopes of the error lines which 
indicate the numerical order of accuracy of the schemes for 


this particular geometry. Of course the grid is not uniform 
in the y direction so the accuracy calculation based on the x 
direction spacing is only approximate. 

5.7. Turbulent Shock Wave Impingement 

The next two-dimensional test case is a shock wave 
impinging on an equilibrium turbulent boundary layer. The 
case is studied experimentally by Reda and Murphy 18 and is 
shown schematically in Fig. 17. A shock generator is 
inclined 13° to the incoming Mach 2.9 freestream. The 
Reynolds number is 5.73xl0 5 per cm, the temperature is 
108° K, and the incoming boundary layer height is 1.694 
cm. The wall is measured to have a temperature of 29 1° K. 
At these conditions, the impinging shock wave creates an 
adverse pressure gradient sufficient to separate the turbulent 
boundary layer. The separation region appears as an 
increase in the thickness of the boundary layer and a shock 
wave is formed off of its leading edge. The flow then 
expands around the separation region and a second shock 
wave is formed at reattachment. Also shown in Fig. 17 is 
the location of the theoretical inviscid shock wave 
impingement location denoted by xj. The results are 
presented with xj as a reference point. This flowfield is 
measured by Reda and Murphy with two different 
experimental setups. In the first set-up, designated by Gl, 
the shock generator spans the width of the channel thereby 
encompassing the sidewall boundary layers. The second 
setup, designated by G2, includes side plates to remove the 
effect of the sidewall boundary layers in order to isolate the 
shock-boundary layer interaction. It is found that including 
side plates reduces the streamwise scale of the interaction 
and lowers the pressure level of incipient separation. The 
flowfields of both setups are inherently three-dimensional in 
nature. It is found that as the shock generator angle is 
increased, separation occurs at the comers of the 
configuration and spreads symmetrically toward the 
centerline. At a shock generator angle of 13°, the flow is 
separated across the channel and a nearly uniform pressure 
distribution is measured. However, the uniformity of the 
pressure distribution across the span is not a sufficient 
condition for two-dimensionality of the flowfield. 
Consequently, comparisons of computed two-dimensional 
results to the experimental data must be made with caution. 

The flowfield is computed using a 101x81 grid. The 
points are equally spaced in the streamwise direction and are 
geometrically stretched away from the wall. The minimum 

wall spacing is near y* = 1. The turbulence is simulated 
using the one-equation turbulence model of Edwards and 
McRae 19 . The first-order grid aligned and CCNO pressure 
contours are shown in Figs. 18 and 19, respectively. The 
grid aligned pressure contours show a smeared incident and 
reflected shock wave. The interaction region does not 
indicate separation of the boundary layer. The CCNO 
pressure contours show the shock waves to be captured with 
less smearing. Furthermore, the pressure contours indicate 
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an interaction region typical of separated flows i.e. both a 
separation and reattachment shock wave are visible. A 
comparison of the skin friction along the wall is plotted in 
Fig. 20. Also shown are the experimentally measured 
separation and reattachment locations. The skin friction 
results verify that the gnd aligned solution does not predict 
a separation and that the CCNO solution does predict a 
separation. Despite the CCNO scheme's prediction of the 
separation zone, both solutions produce similar pressure 
distributions that are not in favorable agreement with the 
experimental data as shown in Fig. 21. The initiation of 
the pressures rise is not predicted well and the shape of the 
distribution is only in qualitative agreement. Although the 
pressure distributions of the two first-order solutions are 
very similar, clearly the CCNO scheme more accurately 
predicts the character of the flowfield. The shock waves are 
captured with less smearing and the separation zone is 
reproduced. 

The grid aligned solution is run for 5,000 iterations at a 
CFL of 5,000 consuming 4.04 minutes of CPU on the 
Cray Y-MP/C90. The CCNO solution is computed with 
the rotation directions based on the pressure gradient. The 
rotation angles are damped to zero below the sonic line. It 
is found that rotating the stencil below the sonic line in the 
separation zone severely inhibits convergence. The solution 
is first computed in grid aligned mode for 1,000 iterations at 
a CFL of 5,000 to allow the incident shock wave to set up. 
The rotation angles are then computed at a CFL of 50 up to 
iteration 9,000 upon which the angles are frozen. The CFL 
is then reduced to 5 and the solution is computed to 
iteration 16,000. The total CPU time is 15.73 minutes on 
the Cray Y-MP/C90. 

The shock impingement case is computed to second- 
order accuracy. The grid aligned pressure contours are 
shown in Fig. 22 while the CCNO pressure contours are 
shown in Fig. 23. It appears that the CCNO solution is 
capturing the incident and reflected shock waves with greater 
clarity. However, examination of the skin friction 
distribution of Fig. 24 reveals the two solutions predict 
very similar results. Both second-order solutions over 
predict the size of the separation zone. The pressure 
distributions through the interaction zone are shown in Fig. 
25. The two schemes produce similar predictions which are 
in good agreement with the G1 data at the initial pressure 
rise. The post interaction pressure levels are in good 
agreement with the experimental data as well. It is difficult 
to determine which of the two experimental data sets are 
best used for a two-dimensional comparison. It is the 
intention of the experiment to reduce the three- 
dimensionality of the flowfield by cutting off the fully 
developed turbulent sidewall boundary layers with the side 
plates. Therefore, one anticipates that the G2 data is best 
for a two-dimensional comparison. However, adding the 
side plates generates laminar sidewall boundary layers that 
are more prone to separation than their turbulent 
counterparts. Thus, it is unclear which set of data is best to 


compare against the computational results. Moreover, the 
effectiveness of the turbulence model on this shock- 
impingement case is uncertain. Regardless of the 
comparison to the experimental data, it is important to 
discover that the CCNO solution does not significantly alter 
the flowfield predictions to second-order accuracy. The 
pressure contours suggest that the CCNO algorithm yields a 
qualitatively better solution in the inviscid region of the 
flowfield, yet the wall predictions are quite similar. 

The second-order grid aligned scheme is run for 19,000 
iterations consuming 15.31 minutes of CPU on the Cray 
Y-MP/C90. The CCNO solution is computed with the 
same rotation angle calculation as the first-order solution. 
The computation is run for 16,500 iterations in 19.07 
minutes. The rotation angles are frozen at iteration 10,000 
and are shown in Fig. 26. 

5.3 Inviscid. Mach 2.8. Channel Flow 

The first three-dimensional test case is Mach 2.8 
inviscid flow in a four wall channel containing a ramp. 
This test case is taken from Rumsey 6 and the geometry is 
shown in Fig. 27. The foot of the ramp is at an oblique 
angle to the oncoming freestream and the inclination of the 
ramp varies across the domain. The side boundaries are all 
solid walls. With an inflow Mach number of 2.8, a truly 
three-dimensional shock surface develops within the channel 
thus providing a challenging test of the three-dimensional 
rotated scheme. 

The case is run on a coarse grid with 41 points in the 
stream wise direction and crossflow planes of 17x17. All 
grid points are evenly spaced along boundary line segments. 
The interior points are linearly interpolated between 
opposite faces. One selected plane of the grid in each 
coordinate direction is shown in Fig. 28. The solution is 
computed to first-order accuracy with the grid aligned 
algorithm. The resulting Mach number contours on three 
selected planes interior to the domain are shown in Fig. 29. 
The contours for this and all subsequent contour plots range 
from 1.4 to 2.8 in increments of 0.05. The predominant 
feature is the ramp shock wave that reflects off of the top 
wall and exits the domain as shown in the j plane view. 
The k and i planes show variations across the domain 
demonstrating the three-dimensional nature of the solution. 
However, the first-order grid aligned results are highly 
smeared and secondary flow features are difficult to discern. 

The solution is computed using the three-dimensional 
CCNO algorithm with the pressure gradient selected as the 
preferred direction. The resulting Mach contours are shown 
in Fig. 30. Clearly, the shock system is captured with less 
smearing using the CCNO algorithm. Not only is the 
primary shock resolved to a higher degree, but secondary 
shock structures that are smeared out with the grid aligned 
scheme are now evident. For example, the k plane view in 
Fig. 30 shows a secondary shock wave emanating from a 
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triple point midway across the ramp shock wave. This 
feature is completely smeared out in the grid aligned 
solution. 

The first-order results of the CCNO algorithm are 
encouraging. Recall that few multi-dimensional upwind 
schemes have been developed and tested in three dimensions. 
Although certainly it can be predicted that aligning the Roe 
solver across a planar shock surface in three dimensions will 
reduce dissipation, it is not certain that any improvement at 
all can be expected by aligning a single coordinate axis in a 
preferred direction when crossing oblique shock waves exist 
such as in this test case. Thus, the first-order results 
present promise for the three-dimensional algorithm. 

In order to quantify the error, a fine grid solution is 
computed and considered as the "exact" solution of an error 
norm calculation. The fine grid mesh spacing is one fourth 
that of the 41x17x17 grid thus yielding a 161x65x65 grid. 
The solution is computed to second-order accuracy with the 
grid aligned scheme. The resulting Mach contours are 
shown in Fig. 31. Evident in the fine grid solution is the 
curvature of the shock waves in the i and k planes as a 
result of the increase in shock strength created by the 
increase in ramp inclination across the domain. Moreover, 
the shock wave reflection off of the side walls creates a 
triple point evident in both the i and k views. A 
comparison of the first-order solutions to the fine grid 
solution is provided in Fig. 32 where a pressure survey 
through the center of the domain in the streamwise direction 
is shown. The pressure survey demonstrates the improved 
shock capturing ability of the CCNO algorithm. 

The L 2 -norm of the mass equation of the first-order grid 
aligned solution based on the fine grid solution is computed 
to be 1.889xl0' 3 . The first-order CCNO error is calculated 
to be 1.1314xl0' 3 . The ratio of the grid aligned to rotated 
error is 1.670. Recall that in the two-dimensional in viscid 
channel flow problem the error ratio is nearly 2. This is an 
indication that the improvement between grid aligned and 
rotated schemes is greater in two dimensions than in three 
dimensions to first-order accuracy. 

The solution on the coarse grid is computed to second- 
order accuracy. The grid aligned Mach contours are shown 
in Fig. 33 while the CCNO Mach contours are shown in 
Fig. 34. From comparison between these two figures, it is 
seen that the rotated scheme improves the three-dimensional 
shock capturing. The CCNO result shows the triple point 
in the i and k planes while the grid aligned scheme only 
hints at its existence. Furthermore, a pressure survey 
through the middle of the domain, as shown in Fig. 35, 
also indicates improved accuracy. The mass equation error 
norm (based on the fine grid solution) of the second-order 
grid aligned scheme is calculated to be 9.566X10" 4 while the 
CCNO error is 7.339X10 -4 . The ratio of the error norms is 
1.30. The improvements of the second-order accurate three- 
dimensional scheme are not as marked as the first-order 


results. However, this is still an encouraging result. The 
multi-dimensional wave model of Rumsey 6 showed no 
improvements to second-order accuracy. Furthermore, keep 
in mind that only a single preferred direction is defined and 
that the rotation is complete when a coordinate axis is 
aligned in that direction. There exists an additional degree 
of freedom in the plane normal to the preferred direction 
where an additional rotation can be defined. Then one of the 
two remaining coordinate axes could be aligned in some 
secondary direction. In this current formulation, this option 
is not explored and thus additional improvements may be 
possible. 

The convergence histories of the coarse grid solutions 
are presented in Fig. 36. The first-order grid aligned 
solution is run at a CFL of 10,000 and reaches a 
convergence level of 10'^ in 170 iterations and 11.4 
seconds on the Cray Y-MP/C90. The second-order grid 
aligned solution is also run with a CFL of 10,000 and 
reaches convergence in 369 iterations taking 23.8 seconds. 
The CCNO solutions are both run at a CFL of 10 with the 
preferred direction aligned with the pressure gradient. For 
the first-order CCNO solution the angles are frozen when 
the residual reaches 9.x 1 O' 4 . Convergence is reached in 297 
iterations and 41 seconds of CPU. The second-order CCNO 
rotation angles are frozen after 150 iterations. Convergence 
is achieved in 593 iterations and 107 seconds. 

Consistent with the two-dimensional inviscid results, it 
is seen that the CCNO algorithm takes longer to converge 
than the grid aligned scheme. This is a result of the 
decrease in numerical dissipation and a poor temporal 
linearizaton of the rotated fluxes. Recall that the implicit 
side does not account for the rotation of the fluxes. It is 
also seen that the CCNO algorithm is significantly more 
computationally expensive than the grid aligned scheme. 
The grid aligned scheme (both first- and second-order) are 
running at approximately 6.3 micro-seconds per point per 
iteration. The first- and second-order CCNO algorithms are 
running at 13.5 and 17.6 micro-seconds per point per 
iteration, respectively. A large portion of the increase is 
due to the second-order interpolation because the 
computation of the interpolation indices does not optimize 
well on the Cray Y-MP/C90. The other intensive 
computation is the three-dimensional transformation of the 
rotated fluxes. An option to reduce this computation is to 
store the results of the transformation after the angles have 
been frozen. The storage requirements are an additional 18 
real numbers per point. Furthermore, the indices of the 
second-order interpolation could also be stored and this could 
perhaps increase the amount of optimization. Although 
there is a penalty in the required memory, this option would 
significandy reduce the CPU usage since a large portion of 
the calculation is performed with the rotation angles frozen. 
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S 4 T .aminar Intersec ting Wedge Comer Flow 

The three-dimensional CCNO algorithm is now applied 
to a viscous flow simulation. The test case is from the 
experiments of West and Korkegi 20 who study the 
interaction region in the comer of two intersecting 9.5 
wedges in a Mach 3 freestream. At a Reynolds number of 
3.9x10^, the flowfield remains laminar. The interaction of 
the two wedge shock waves forms an oblique comer shock 
wave Originating at the intersection of the wedge and 
comer shock waves are an embedded shock wave and a shear 
surface that terminates at the comer. The flowfield is 
discretized using 41 point in the streamwise direction spaced 
along 1.45 units in the freestream direction. The crossflow 
planes are 51x51 with a minimum spacing of 0.0005 units 
at the wall and then stretched to the outer boundary located 
at 1.75 units. The freestream temperature is 97.5 Kelvin 
and the wall temperature is held constant at 2.8 times the 
freestream temperature. 

The solution of the intersecting wedge configuration is 
computed to first-order accuracy using both the grid aligned 
and CCNO algorithms. The rotation orientation of the 
CCNO algorithm is based on the Mach number gradient. 
Pitot pressure contours on three crossflow planes of the grid 
aligned and CCNO solutions are shown in Figs. 37 and 38, 
respectively. The first-order grid aligned solution shows the 
two wedge shock waves to be joined by a curved shock 
wave in the comer. The two embedded shock waves are 
seen to be smeared. The first-order CCNO contours show 
the wedge shock waves to be captured with less smearing. 
Furthermore, the union of the two wedge shock waves 
occurs through a comer shock wave and not a continuous 
transition as predicted by the grid aligned scheme. Also 
evident in the CCNO solution is the slip surfaces 
terminating at the comer. These slip surfaces are not seen 
in the grid aligned solution at the current contour settings. 

As in the two-dimensional calculations, the inviscid 
portion of the three-dimensional flowfield is captured with 
greater clarity with the CCNO scheme to first-order 
accuracy. The previous pitot pressure contours reveal a 
qualitative improvement of the CCNO scheme as compared 
to the grid aligned. A quantitative comparison of the two 
solutions is contained in Fig. 39 where the wall pressure 
distribution in the transverse direction is compared to the 
experimental data. The ordinate of the plot is 
nondimensionalized by the x location of the survey i.e. x = 
1.45. Also, the abscissa reference value, p is 1.98 
times p„. Figure 39 shows that both first-order solutions 
are in general agreement with the experimental data. Both 
solutions are seen to under predict the pressure in the comer 
inside of the imbedded shock wave. This is a result of the 
coarseness of the grid. Prior to the pressure drop across the 
embedded shock wave , the CCNO solution is shown to 
slightly indicate the pressure rise between z/x = 0.2 and z/x 
= 0.3 whereas the grid aligned solution predicts a level 


pressure plateau. Both solutions are shown to under predict 
the position of the embedded shock wave along the z/x axis. 
Across the embedded shock wave, the smearing of the grid 
aligned solution is evident. Outboard of the embedded 
shock wave, the two solutions show a similar prediction 
which is in general agreement with the experimental data. 
Overall, the prediction of the wall pressure distributions is 
shown to be only modestly improved with the CCNO 
algorithm despite the qualitative improvements in the 
flowfield. Recall that the first-order two-dimensional results 
show that the improvements in the inviscid portion of the 
flowfield yielded more accurate wall predictions. In this 
three-dimensional test case, this improvement is not 
realized. 

The solution is also computed to second-order accuracy. 
The resulting pitot pressure contours of the grid aligned and 
rotated algorithms are shown in Figs. 40 and 41, 
respectively. Comparing the two second-order solutions, 
only small differences in the character of the flowfield are 
detected. The CCNO solution is shown to capture the 
embedded shock waves with slightly more clarity. This is 
best seen in the middle of the three crossflow planes. 
Comparison of the wall pressure distribution is contained in 
Fig. 42. Evident in the comer region are low amplitude 
pressure oscillations of the CCNO solution. These 
oscillations are a result of the rotation in the boundary 
layer. As mentioned in the discussion of the two- 
dimensional turbulent shock impingement case, the 
solution is best computed by damping out the rotation 
below the sonic line. In this calculation, the angles are not 
damped and the pressure oscillations develop. The largest 
difference in the wall pressure distributions is seen 
immediately outboard of the embedded shock wave. The 
grid aligned algorithm is shown to under predict the pressure 
in the crossflow reattachment region. The under prediction 
of the pressure through this region is characteristic of 
upwind forumulations as can be seen in the results of Rudy, 
et. al. 21 . The CCNO solution is shown to be in better 
agreement with the experimental data in this region. 
Interestingly enough, the study of Rudy, et. al. shows that 
only a central difference formulation (MacCormack’s 
scheme) did not predict the pressure dip suffered by the 
upwind schemes. In the case of the CCNO solution, it 
appears that realignment of the upwind operator in this 
interaction region, reproduces results characteristic of a 
central difference formulation. 

The first-order grid aligned solution reaches convergence 
in 3,707 iterations using 52 minutes of processor time on 
the Cray Y-MP/C90. The second-order grid aligned 
solution is run for 9,000 iterations and 126 minutes of 
computer time reaching a five order of magnitude 
convergence. The CCNO solutions are computed with the 
rotation based on the Mach number gradient. The rotation 
angles are frozen after 3,000 iterations. The first-order 
solution is run for 10,000 iterations consuming —50 
minutes of computer time. The second-order solution is run 
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for 9,000 iterations which uses 273 minutes of computer 
time. The two rotated solutions converge four to five orders 
of magnidue. 

6 0 Conclusion 

This study has presented the development of a rotated 
upwind algorithm for the numerical solution of the Euler 
and Navier-Stokes equations in both two and three 
dimensions. The two-dimensional algorithm is shown to 
predict an inviscid channel flow solution with greater 
accuracy than traditional grid aligned upwinding to both 
first- and second-order accuracy. The improvements to the 
second-order accurate results are shown to be as great as the 
first-order accurate solutions. Futhermore, the accuracy 
gains are maintained as the grid is refined. The calculation 
of a turbulent shock impingement flowfield shows that the 
new algorithm improves the shock wave capturing of the 
inviscid portion of the flowfield to both first- and second- 
order accuracy. However, only the first-order accurate 
scheme shows improved agreement with experimental wall 
pressure and skin friction distributions. 

In three dimensions, a series of rotation matrices is 
developed that aligns a single coordinate axis in a preferred 
direction. Furthermore, the rotation matrices are defined 
such that the orientation of the rotated coordinate axes with 
respect to the grid axes is always known. The three- 
dimensional algorithm is applied to an inviscid shock 
reflection problem. The new algorithm is shown to predict 
the flowfield with greater accuracy to both first- and second- 
order accuracy. However, the improvements to the three- 
dimensional solution are not as great as those occurring in 
the two-dimensional inviscid solution. The new algorithm 
is shown to be sufficiently robust to compute the flow in 
the comer of intersecting wedges. To both first- and second- 
order accuracy, the inviscid features of the flowfield are 
shown to be qualitatively improved with the rotated 
algorithm. However, surface pressure predictions are only 
marginally improved. 

Both the two- and three-dimensional algorithms are 
shown to consume 2-4 times as much CPU time as the grid 
aligned solutions. Further possibilities exist for improving 
the current algorithm. A better temporal linearization of the 
rotated fluxes may improve the convergence rate. Also, 
simplified interpolation schemes may produce similar 
results at less cost. 
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b. Physical Space 

Figure L) Cell-Center Non-Orthogonal Flux Strategy. 



Cell Center 

First Order Line 
of Interpolation 


Figure 2.) Second-Order Interpolation Stencil in the First 
Quadrant of Computational Space for the 
CCNO Scheme. 
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Figure 3.) Orientation of Rotated Coordinate System for 
Type I Rotation. 



Figure 4.) Orientation of Rotated Coordinate S ystem for 
Type II Rotation. 


c 



Figure 5.) First-Order Interpolating Stencil of +£' 
Coordinate Ray under Type II Rotation. 


f 



Scheme. 



Figure 7.) Inviscid, Mach 2, Channel Flow, 61x21 Grid. 



Figure 8.) Inviscid, Mach 2, Channel Flow Mach 

Contours, First-Order Grid Aligned Scheme, 
61x21 Grid, Increment = 0.05. 


Figure 9.) Inviscid, Mach 2, Channel Flow Mach 

Contours, First-Order CCNO Scheme, 61x21 
grid. Rotation based on Vp , Increment = 0.05. 
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Figure 10.) Inviscid, Mach 2, Channel Flow Steady State 
Upwind Directions, First-Order CCNO Scheme, 
61x21 grid. Rotation based on V 7 /? . 



Figure 14.) Inviscid, Mach 2, Channel Flow Mach 

Contours, Second-Order CCNO Scheme, 61x21 
Grid, Rotation based on V 7 /? , Increment=0.05. 




Figure 11.) Inviscid, Mach 2, Channel Flow Pressure 

Survey along y=10 computational line, First- 
Order Comparison, 61x21 Grid. 


Figure 15.) Inviscid, Mach 2, Channel Flow Pressure 

Survey along y= 10 computational line, Second- 
Order Comparison, 61x21 Grid. 



Figure 12.) Inviscid, Mach 2, Channel Flow Mach 

Contours, Second-Order Grid Aligned Scheme, 
385x129 Grid, Increment=0.05. 



Figure 13.) Inviscid, Mach 2, Channel Flow Mach 

Contours, Second-Order Grid Aligned Scheme, 
61x21 Grid, Increment = 0.05. 
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Figure 16.) Inviscid, Mach 2, Channel Flow Error 
Comparison. 



Figure 17.) Turbulent Shock Wave Impingement Geometry 
and Flowfield Schematic. 
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Figure 18.) Turbulent Shock Wave Impingement Pressure 
Contours, First-Order Grid Aligned Scheme. 



Figure 22.) Turbulent Shock Wave Impingement Pressure 
Contours, Second-Order Grid Aligned Scheme. 



Figure 19.) Turbulent Shock Wave Impingement Pressure 
Contours, First-Order CCNO Scheme. 



Figure 23.) Turbulent Shock Wave Impingement Pressure 
Contours, Second-Order CCNO Scheme. 




Figure 20.) Turbulent Shock Wave Impingement Wall Skin 
Friction Distribution, First-Order Comparison. 



Figure 21.) Turbulent Shock Wave Impingement Wall 

Pressure Distribution, First-Order Comparison. 


Figure 24.) Turbulent Shock Wave Impingement Wall Skin 
Friction Distribution, Second-Order 
Comparison. 



Figure 25.) Turbulent Shock Wave Impingement Wall 
Pressure Distribution, Second-Order 
Comparison. 
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Figure 26,) Turbulent Shock Wave Impingement Steady 

State Upwind Directions, Second-Order CCNO 
Scheme, 



Figure 27.) Inviscid, Mach 2.8, Channel Row Geometry. 


a. k = 9 plane. b. i = 21 plane. 



c. j = 9 plane. 


Figure 28.) Inviscid, Mach 2.8, Channel Row, Selected 
Grid Planes, 41x17x17 Grid. 



a. k = 8 plane. b. i = 20 plane. 



c. j = 8 plane. 


Figure 29.) Inviscid, Mach 2.8, Channel Row Mach 

Contours, First-Order Grid Aligned Scheme, 
41x17x17 Grid. 



a. £ = 8 plane. b. i = 20 plane. 



c. ; = 8 plane. 

Figure 30.) Inviscid, Mach 2.8, Channel Row Mach 
Contours, First-Order CCNO Scheme, 
41x17x17 Grid. 
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a. k = 32 plane. b. i = 80 plane. 



c. j = 32 plane. 

Figure 31.) Inviscid, Mach 2.8, Channel Flow Mach 

Contours, Second-Order Grid Aligned Scheme, 
161x65x65 Grid. 




c. j = 8 plane. 

Figure 34.) Inviscid, Mach 2.8, Channel How Mach 
Contours, Second-Order CCNO Scheme, 
41x17x17 Grid. 



Figure 32.) Inviscid, Mach 2.8, Channel Flow Pressure 

Survery along j = 8, k = 8 computational line, 
First-Order Comparison, 41x17x17 Grid. 



Figure 35.) Inviscid, Mach 2.8, Channel Flow Pressure 

Survery along j = 8, k = 8 computational line, 
Second-Order Comparison, 41x17x17 Grid. 



a. k = 8. 


b. 



i = 20 plane. 



c. j = 8 plane. 


Figure 33.) Inviscid, Mach 2.8, Channel Flow Mach 

Contours, Second-Order Grid Aligned Scheme, 
41x17x17 Grid. 



Figure 36.) Inviscid, Mach 2.8, Channel How Convergence 
History, 41x17x17 Grid. 
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Figure 37 ) Laminar Intersecting Wedge Comer Flow Pitot Figure 40.) Laminar Intersecting W edge Comer Flow Pitot 
Pressure Contours, First-Order Grid Aligned Pressure Contours, Second-Order Grid Aligned 

Scheme. Scheme. 



Figure 38.) Laminar Intersecting Wedge Comer Flow Pitot Figure 41.) Laminar Intersecting Wedge Comer Flow Pitot 
Pressure Contours, First-Order CCNO Scheme. Pressure Contours, Second-Order CCNO 

Scheme. 
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Figure 39.) Laminar Intersecting Wedge Comer Flow Wall 
Pressure Distribution, First-Order Comparison. 
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Figure 42.) Laminar Intersecting Wedge Comer Flow Wall 
Pressure Distribution, Second-Order 
Comparison. 
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